[Fluidentity13]

Need Help !! How do you solve linear System of Equations(Solver), with a Singular or Ill-Conditioned Matrix?

The linear system of equation is based on Pressure corrector Step 2 from incompressible SIMPLE algorithm, & I would be implementing same for converging loops of SIMPLE for X & Y velocity.
I am stuck on the ill-conditioned matrix (cond = 1e16, please only consider first one in the image) & can't seem to find a solution to this problem. Can anybody point me in the correct direction. I've picked up the derivation from Roy-Chowdhury/Patankar, on a structured 10/10 Mesh.


Can a linear solver(Python based) do this, or should I go about coding steepest descent or conjugate gradient, as I did with a relaxation sample problem before.

I've attached the images of the A matrix, cond number range of matrix, & the warning message I receive in the console during run.

Also here is a code snippet of the Pressure Correction Matrix formulation, if anyone's interested.

App =np.zeros((ny, nx))
Ape =np.zeros((ny, nx))
Apw =np.zeros((ny, nx))
Apn =np.zeros((ny, nx))
Aps =np.zeros((ny, nx))
Bpp =np.zeros((ny, nx))
for j in range (0, nx):
for i in range (0, ny):
if(j<nx-1):
Ape[i, j] = rho*dy**2/Aup[i, j]
else:
Ape[i, j] = 0
if(j>0):
Apw[i, j] = rho*dy**2/Aup[i, j-1]
else:
Apw[i, j] = 0
if(i>0):
Apn[i, j] = rho*dx**2/Avp[i-1, j]
else:
Apn[i, j] = 0
if(i<ny-1):
Aps[i, j] = rho*dx**2/Avp[i, j]
else:
Aps[i, j] = 0
Bpp[i, j] = rho*dy*(uxs[i, j]-uxs[i, j+1]) + rho*dx*(uys[i+1, j]-uys[i, j])
App[i, j] = Ape[i, j] + Apw[i, j] + Apn[i, j] + Aps[i, j]

App_ = np.flip(App, 0)
Ape_ = np.flip(Ape, 0)
Apw_ = np.flip(Apw, 0)
Apn_ = np.flip(Apn, 0)
Aps_ = np.flip(Aps, 0)
Bpp_ = np.flip(Bpp, 0)
pc_ = np.zeros((ny, nx))

pcm_ = pc_.flatten()
Bppm_ = Bpp_.flatten()

A = np.zeros((ny*nx, ny*nx))
c=0
for i in range (0, ny):
for j in range (0, nx):
A[c, c] = App_[i, j]
if ((c+1)%nx!=0):
A[c, c+1] = -Ape_[i, j]
if(c%nx!=0 or c!=0):
A[c, c-1] = -Apw_[i, j]
if (c<(ny-1)*nx):
A[c, c+nx] = -Apn_[i, j]
if(c>nx-1):
A[c, c-nx] = -Aps_[i, j]
c+=1

Please advice on posting with proper code Indentation too. I'll correct this.

[FMDenaro]

Let me generalize the issue. The fact that the pressure matrix is singular is a consequence of the Neumann BCs and expresses the fact that the pressure solution is determined up to a function of time.
A solution is determined even using an iterative method, provided that the BCs fulfill the compatibility relation.

I suppose that a non-convergence is a signal of wrong BCs

[Fluidentity13]

So, you mean the formulation of the pressure correction matrix must be wrong, and it's not commonplace for illconditioned matrix, requiring special attention. Sorry, I am very new to this. And rechecking so many times, made me believe something else needed to be done for such a matrix.

Coming on to the boundary conditions. I am solving this in a staggerred system. And the pressure corrector system contains as many equations as there are cells in the mesh. 10 X 10. Owing to the staggerred scheme I am able to include both side u_star terms in the rhs B matrix. Can u take a look at the Matrix Formulation image I've attached for reference. I am not sure what to do next!!

[FMDenaro]

Originally Posted by Fluidentity13

So, you mean the formulation of the pressure correction matrix must be wrong, and it's not commonplace for illconditioned matrix, requiring special attention. Sorry, I am very new to this. And rechecking so many times, made me believe something else needed to be done for such a matrix.

Coming on to the boundary conditions. I am solving this in a staggerred system. And the pressure corrector system contains as many equations as there are cells in the mesh. 10 X 10. Owing to the staggerred scheme I am able to include both side u_star terms in the rhs B matrix. Can u take a look at the Matrix Formulation image I've attached for reference. I am not sure what to do next!!

How do you prescribe the BCs for the pressure?

[Fluidentity13]

I haven't, prescribed anything. Initial conditions I've provided. I'll have to revisit the chapters once again, I guess. Is the boundary condition treated by equating the last row to a fixed value, or to the second last row.

Sorry I got really confused, devising something for the staggerred scheme, but I guess it's getting better now with your help.

[FMDenaro]

Originally Posted by Fluidentity13

I haven't, prescribed anything. Initial conditions I've provided. I'll have to revisit the chapters once again, I guess. Is the boundary condition treated by equating the last row to a fixed value, or to the second last row.

Sorry I got really confused, devising something for the staggerred scheme, but I guess it's getting better now with your help.

have you already searched in the forum for similar post?

[Fluidentity13]

No, I guess I need information on specifying boundary conditions.

Only one wall needs Pressure boundary right? And the whole matrix would resize because of that? Also for X-Velocity, I have provided Left and Right Boundary Cond, and for Y Velocity, Top And Bottom boundary are equated to Zero. Does this all sound good to you?

[FMDenaro]

Originally Posted by Fluidentity13

No, I guess I need information on specifying boundary conditions.

Only one wall needs Pressure boundary right? And the whole matrix would resize because of that? Also for X-Velocity, I have provided Left and Right Boundary Cond, and for Y Velocity, Top And Bottom boundary are equated to Zero. Does this all sound good to you?

Poisson equation,Neumann BCs

[Fluidentity13]

Thank a lot for diverting me to that enlightening thread. I picked up that only applying Nuemann Boundary condition can tease out a solution, & wanted to give it a try, assigning the right boundary that condition. Instead of adding nuemann terms to the source term, following changes were made to the A matrix in Ax=b:

1. Bottom row of each submatrix removed which contained rightmost cell coefficients. (Similar changes to p & B matrix.
2. New bottom row of each submatrix changed to App-Ape from plain App. (Neumann terms)

Anyways, that did not work for me still giving that error for inaccurate solution, cond (1e17). Does my implementation look good to you, must be a problem in the code?

Assigning the top-left most cell dirichlet boundary did the job for the con number (1040). There's still huge divergence in the whole algorithm somewhere though, I am sure there's some mistake somewhere.

To summarise, the mistake I did was not providing any boundary condition to the pressure corrector equation whatsoever. It's recommended to fix a cell (mostly to 0), & update the Ax=b matrix based on that.

[FMDenaro]

Originally Posted by Fluidentity13

Thank a lot for diverting me to that enlightening thread. I picked up that only applying Nuemann Boundary condition can tease out a solution, & wanted to give it a try, assigning the right boundary that condition. Instead of adding nuemann terms to the source term, following changes were made to the A matrix in Ax=b:

1. Bottom row of each submatrix removed which contained rightmost cell coefficients. (Similar changes to p & B matrix.
2. New bottom row of each submatrix changed to App-Ape from plain App. (Neumann terms)

Anyways, that did not work for me still giving that error for inaccurate solution, cond (1e17). Does my implementation look good to you, must be a problem in the code?

Assigning the top-left most cell dirichlet boundary did the job for the con number (1040). There's still huge divergence in the whole algorithm somewhere though, I am sure there's some mistake somewhere.

To summarise, the mistake I did was not providing any boundary condition to the pressure corrector equation whatsoever. It's recommended to fix a cell (mostly to 0), & update the Ax=b matrix based on that.

Fixing a value is a way to get one solution because you fixed the constant value. But that does not mean you have the correct solution.

Boundary conditions are required for the mathematical pressure problem to be wel posed. You need then to discretize those BCs accordingly.

[vyas.manutd89]

I am facing a similar issue. I am simulating a steady state 2D channel flow in my code with velocity inlet and pressure outlet B.C`s. For the pressure correction equation I have imposed Neumann B.C at the top and bottom and zero correction at the inlet & outlet. What i am seeing is that the sparse matrix for pressure correction is becoming ill-conditioned (almost singular) due to which I am getting extremely high values for pressure correction with the 1st iteration itself.



Would be very interested to know how to solved the same issue that.

[Fluidentity13]

Zero Correction at inlet ? Shouldn't be that way for channel flow, if I understand your inquiry right. Only for Outlet you've to fix the Pressure value.

Also, staggered or collocated grid.

You can reference my code if you want.



2D_Panel-CFD v0.0.1 Release

[vyas.manutd89]

I am a little confused with the pressure correction boundary condition for the inlet/inflow condition. Patankar says with inflow there is not need to prescribe a pressure correction B.C.
I was able to obtain the solution using collocated grid & momentum interpolation approach, and i used a Neumann condition at the inlet and walls.


But when trying to replicate in a staggered grid arrangement I am running into problems with a ill-conditioned sparse matrix for solving pressure correction. I guess pressure correction being elliptic, it is important to prescribe the correct B.C`s on all the walls of the domain.

[FMDenaro]

Originally Posted by vyas.manutd89

I am a little confused with the pressure correction boundary condition for the inlet/inflow condition. Patankar says with inflow there is not need to prescribe a pressure correction B.C.
I was able to obtain the solution using collocated grid & momentum interpolation approach, and i used a Neumann condition at the inlet and walls.


But when trying to replicate in a staggered grid arrangement I am running into problems with a ill-conditioned sparse matrix for solving pressure correction. I guess pressure correction being elliptic, it is important to prescribe the correct B.C`s on all the walls of the domain.

In inflow you prescribe the Dirichlet bc for the velocity but the pressure equation requires to prescribe bcs on the whole frontier of the domain. In your case you have non-homogenoeus Neumann BC to prescribe at inflow and walls.