CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Help with Pressure Correction Matrix - Singular - Ill-conditioned Matrix.

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By FMDenaro

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 1, 2021, 07:30
Exclamation Help with Pressure Correction Matrix - Singular - Ill-conditioned Matrix.
  #1
New Member
 
Aman Sangal
Join Date: Mar 2021
Posts: 21
Blog Entries: 2
Rep Power: 3
Fluidentity13 is on a distinguished road
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.



Quote:
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.
Attached Images
File Type: jpg Coefficient matrix for pressure correction.jpg (114.1 KB, 14 views)
File Type: png ill-conditioned matrix .png (45.2 KB, 10 views)
File Type: png Warning Message .png (76.8 KB, 7 views)
File Type: jpg Matrix Formulation.jpg (154.3 KB, 13 views)

Last edited by Fluidentity13; November 1, 2021 at 07:35. Reason: Code Indentation
Fluidentity13 is offline   Reply With Quote

Old   November 1, 2021, 08:04
Default
  #2
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,062
Rep Power: 64
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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
FMDenaro is offline   Reply With Quote

Old   November 1, 2021, 08:16
Default
  #3
New Member
 
Aman Sangal
Join Date: Mar 2021
Posts: 21
Blog Entries: 2
Rep Power: 3
Fluidentity13 is on a distinguished road
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!!
Fluidentity13 is offline   Reply With Quote

Old   November 1, 2021, 08:18
Default
  #4
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,062
Rep Power: 64
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Fluidentity13 View Post
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?
FMDenaro is offline   Reply With Quote

Old   November 1, 2021, 08:24
Default
  #5
New Member
 
Aman Sangal
Join Date: Mar 2021
Posts: 21
Blog Entries: 2
Rep Power: 3
Fluidentity13 is on a distinguished road
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.
Fluidentity13 is offline   Reply With Quote

Old   November 1, 2021, 08:29
Default
  #6
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,062
Rep Power: 64
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Fluidentity13 View Post
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?
FMDenaro is offline   Reply With Quote

Old   November 1, 2021, 08:38
Default
  #7
New Member
 
Aman Sangal
Join Date: Mar 2021
Posts: 21
Blog Entries: 2
Rep Power: 3
Fluidentity13 is on a distinguished road
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?
Fluidentity13 is offline   Reply With Quote

Old   November 1, 2021, 08:40
Default
  #8
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,062
Rep Power: 64
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Fluidentity13 View Post
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 likes this.
FMDenaro is offline   Reply With Quote

Old   November 2, 2021, 03:32
Default
  #9
New Member
 
Aman Sangal
Join Date: Mar 2021
Posts: 21
Blog Entries: 2
Rep Power: 3
Fluidentity13 is on a distinguished road
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.
Fluidentity13 is offline   Reply With Quote

Old   November 2, 2021, 04:12
Default
  #10
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,062
Rep Power: 64
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Fluidentity13 View Post
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.
FMDenaro is offline   Reply With Quote

Old   May 3, 2022, 11:06
Default
  #11
New Member
 
Vyas Duggirala
Join Date: May 2014
Posts: 4
Rep Power: 10
vyas.manutd89 is on a distinguished road
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.
vyas.manutd89 is offline   Reply With Quote

Old   May 3, 2022, 14:07
Default
  #12
New Member
 
Aman Sangal
Join Date: Mar 2021
Posts: 21
Blog Entries: 2
Rep Power: 3
Fluidentity13 is on a distinguished road
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
Fluidentity13 is offline   Reply With Quote

Old   May 4, 2022, 01:46
Default
  #13
New Member
 
Vyas Duggirala
Join Date: May 2014
Posts: 4
Rep Power: 10
vyas.manutd89 is on a distinguished road
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.
vyas.manutd89 is offline   Reply With Quote

Old   May 4, 2022, 13:43
Default
  #14
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,062
Rep Power: 64
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by vyas.manutd89 View Post
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.
FMDenaro is offline   Reply With Quote

Reply

Tags
ill-conditioned matrix, implicit solvers, pressure correction, simple, singular matrix

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
What is difference between static pressure and gauge pressure? aja1345 FLUENT 1 July 20, 2018 20:05
question regarding LES of pipe flow - pimpleFoam Dan1788 OpenFOAM Running, Solving & CFD 37 December 26, 2017 14:42
SIMPLE algorithm does not converge when using old pressure (correction) values andreasp Main CFD Forum 3 February 9, 2016 21:18
Pressure correction equation,Ferziger and Peric Hooman Main CFD Forum 1 July 20, 2010 05:21
Hydrostatic pressure in 2-phase flow modeling (long) DS & HB Main CFD Forum 0 January 8, 2000 15:00


All times are GMT -4. The time now is 11:43.