[AeroSonic]

Dear All,

I have a question regards my Solver, I have build a FORTRAN code to solve the incompressible viscous N--S equations.

at first I used dgetrf and dgetri to get the inverse of the pressure matrix and solve the pressure system by matrix multiplication.

all went well and good, after that I tried to include the PARDISO solver so I can save memory when I use a very large grid. so I am solving the pressure system by PARDISO instead of dgetrf and dgetri and matrix multiplication.


I compared the results at the first few time steps and it was identical, and tried that too at time step =100,1000,5000 and the solution was identical too. the problem appears when the time steps increases, it is like there is a cumulative error which increase with time.

my grid is 20Dx10D , Re=200, Epsilon = 0.05, dt=0.005

there is a sample of the results in attachments.
I can see the difference in numbers grows with time steps and I dont know why, although the solution was identical at the first time steps so there is no mathematical problem in the code.

[AeroSonic]

Boundary condition

inlet (u=1,v=0)
upper and lower ( v=0)
exit ( P at the midline point = 0 )

Initial Coniditons : from inlet

[FMDenaro]

Difficult to say what the problem could be.
I would address the fact you wrote "the inverse of the pressure matrix". In principle, the pressure matrix is singular, how do you treat this issue?

[sbaffini]

Seems like a good case for doing a verification as opposed to a validation:

1) pick up an unsteady case with analytical solution (stokes, taylor vortex, etc.)
2) solve it with different grids and constant CFL, maybe at multiple CFL numbers
3) For each CFL and each grid-time step combination use both solvers
4) Compare both results with the analytical ones

Nothing can beat this at spotting problems and if differences are always below the threshold of the numerical error you can maybe try again with a grid-time step-CFL combination that will bring it below the differences you see. At this point, either the differences will go further down or they won't.

But, beside this, I would maybe first try the following:

1) set up a test with a matrix A like the one from your pressure equation. If you work in 2D, parameterize it with both dx and dy and make different attempts with different values
2) assign yourself a solution x (say, all ones) and compute the resulting rhs = Ax
3) solve the system with both solvers using A and rhs
4) compare results x and try to find out what's wrong and where
5) If it turns out that it is not your fault, just pick the one solver that actually works for you and maybe submit the test to the developers of the one that doesn't

[AeroSonic]

Originally Posted by FMDenaro

Difficult to say what the problem could be.
I would address the fact you wrote "the inverse of the pressure matrix". In principle, the pressure matrix is singular, how do you treat this issue?

I put the BC before getting the Inverse, I get the inverse once outside the time loop.

[FMDenaro]

Originally Posted by AeroSonic

I put the BC before getting the Inverse, I get the inverse once outside the time loop.

What kind of BC are you prescribing?
Then, save a small matrix like a 4x4 and check the inverse

[AeroSonic]

Originally Posted by sbaffini

Seems like a good case for doing a verification as opposed to a validation:

1) pick up an unsteady case with analytical solution (stokes, taylor vortex, etc.)
2) solve it with different grids and constant CFL, maybe at multiple CFL numbers
3) For each CFL and each grid-time step combination use both solvers
4) Compare both results with the analytical ones

Nothing can beat this at spotting problems and if differences are always below the threshold of the numerical error you can maybe try again with a grid-time step-CFL combination that will bring it below the differences you see. At this point, either the differences will go further down or they won't.

But, beside this, I would maybe first try the following:

1) set up a test with a matrix A like the one from your pressure equation. If you work in 2D, parameterize it with both dx and dy and make different attempts with different values
2) assign yourself a solution x (say, all ones) and compute the resulting rhs = Ax
3) solve the system with both solvers using A and rhs
4) compare results x and try to find out what's wrong and where
5) If it turns out that it is not your fault, just pick the one solver that actually works for you and maybe submit the test to the developers of the one that doesn't

I believe the code is running correctly at the first few time steps, as I compared the result at the first time step and after 10 time steps and after 100 time steps, the pressure and velocity results are identical, the first time step is very very similar.

attached for example the result at specific point after 2 time steps

[AeroSonic]

Originally Posted by FMDenaro

What kind of BC are you prescribing?
Then, save a small matrix like a 4x4 and check the inverse

Dirichlet BC, P=0 at the mid point of the exit section.
I am solving the problem in 2D.

[FMDenaro]

Originally Posted by AeroSonic

Dirichlet BC, P=0 at the mid point of the exit section.
I am solving the problem in 2D.

Ok, with Dirichlet BC the matrix is not singular.
Now I suggest to check for the divergence-free constraint in each cell.
It is satisfied at machine precision everywhere?

[AeroSonic]

Originally Posted by FMDenaro

Ok, with Dirichlet BC the matrix is not singular.
Now I suggest to check for the divergence-free constraint in each cell.
It is satisfied at machine precision everywhere?

I am using Finite Element ( Galerkin Least Squares ) so I had to use a stabilization parameter in the continuity as attached

for the second question, how do I know ?

[AeroSonic]

I will solve the lid driven cavity problem and check what is happening, may bw something appear and tell me what is going on