[mihirmakwana6]

Hello Friends!!!

I have written a 2D Navier Stokes solver using Fractional step Method on staggered grid using FVM and tested it for the lid driven cavity problem. I compared my solution with the benchmark results from GHIA ET AL. (1982) JOURNAL OF COMPUTATIONAL PHYSICS VOL. 48, pp.387-411 . I am getting satisfactory results but I am facing an issue with convergence.


My code takes a lot of time to converge when the number of point in x and y direction are more than 22. Hence I decided to use Successive over-relaxation ( SOR ) .

The optimum value for the relaxation factor that I am using is

w_opt = 2 /( 1 + sin(acos(0.5*(cos((22/7)/imax)+cos((22/7)/jmax)))))

where imax,jmax = maximum points in x and y direction

I am using 70% of optimum value in my code.

However with SOR the code takes even more time.

What can be the possible problem with convergence and SOR ?


The equations that I have used are as follows



u_star = u_n + dt * G ( u_i )

where dt = timestep / 2 and G ( u_i ) = G ( u_n ) .... RK2 step 1

dt = timestep and G ( u_i ) = G ( u_(n+1/2) ) .... RK2 step 2

where G = - C + D ( C : convection , D : viscous dissipation )

Within each RK2 step an poisson equation for pressure , P is solved using

Laplacian( P ) = ( density / dt ) Div ( u_star )

where u_n : old time level velocity
u_(n+1/2) : mid time level velocity

Can you please suggest a paper on Fractional step method that uses SOR and also talks about the convergence issue?

Thanks in advance !!

[FMDenaro]

You considered two different issues and the convergence with SOR has nothing to do directly with the FTS method.

The elliptic equation for pressure has a solution (apart a constant) if the compatibility relations is satisfied. Assumed that you checked for this condition, have you fixed a value for the pressure at some location? If so, do not fix any value and let the solver fix its constant.

Then the optimal value you evaluated by the formula how many iterations require? have you checked the number of iterations when you slightly increase and decrease such value? Try do to that always running the first dt and compare. Check also if the divergence of the velocity is satisfied in each cell.

[mihirmakwana6]

Originally Posted by FMDenaro

You considered two different issues and the convergence with SOR has nothing to do directly with the FTS method.

The elliptic equation for pressure has a solution (apart a constant) if the compatibility relations is satisfied. Assumed that you checked for this condition, have you fixed a value for the pressure at some location? If so, do not fix any value and let the solver fix its constant.

Then the optimal value you evaluated by the formula how many iterations require? have you checked the number of iterations when you slightly increase and decrease such value? Try do to that always running the first dt and compare. Check also if the divergence of the velocity is satisfied in each cell.

Thankyou for the reply Filippo Maria Denaro .

which compatibility relation are you talking about? can you please elaborate on it?

I have used a zero gradient at the wall as the Boundary condition for pressure. Before the code runs, I initialized the pressure matrix with the value of zero over the entire domain. Will this initialization create any sort of problem?

Also, If I used the optimal value as given by the formula then the solution blows up. It gives a value of NAN as the solution. I first started with 50% of the optimal value and then based on the amount of time the code takes to reach the convergence ( for steady state as 1e-5), I reached a value of 70% of w_opt as the required value.



For a 32 X 32 grid , Re = 100 , timestep = 0.001s and 70% of w_opt , the code takes

4057 iterations and 44.8 minutes for FOU scheme,
4137 iterations and 54.5 minutes for QUICK scheme and
4104 iterations and 44.5 minutes for CDS scheme.

I run the loop for pressure till the mass source ( DIV(u_star) ) , in each cell is less than 1e-7.

[FMDenaro]

Originally Posted by mihirmakwana6

Thankyou for the reply Filippo Maria Denaro .

which compatibility relation are you talking about? can you please elaborate on it?

I have used a zero gradient at the wall as the Boundary condition for pressure. Before the code run, I initialized the pressure matrix with the value of zero over the entire domain. Will this initialization create any sort of problem?

Also, If I used the optimal value as given by the formula then the solution blows up. It gives a value of NAN as the solution. I first started with 50% of the optimal value and then based on the amount of time the code takes to reach the convergence ( for steady state as 1e-5), I reached a values of 70% of the optimum value.

for a 32 X 32 grid , Re = 100 , timestep = 0.001s and 70% of w_opt , the code takes

4057 iterations and 44.8 minutes for FOU scheme,
4137 iterations and 54.5 minutes for QUICK scheme and
4104 iterations and 44.5 minutes for CDS scheme.

I run the loop for pressure till the mass source ( DIV(u_star) ) , in each cell is less than 1e-7.

I think your FTS code is not properly written...fulfilling the compatibility relation is fundamental to ensure the existence of a solution otherwise the system has no solutions.
After you set homogeneous Neumann boundary conditions, you modified the source term of the elliptic equation?

Starting a guess value = 0 everywhere is not optimal... you can start this way for the first dt, then is better to start from the field at the previous time step.

However, if the iterative solver does not converge for omega = 1 (simple Gauss-Seidel) then you implemented in a wrong way the FTS.

[mihirmakwana6]

Originally Posted by FMDenaro

I think your FTS code is not properly written...fulfilling the compatibility relation is fundamental to ensure the existence of a solution otherwise the system has no solutions.
After you set homogeneous Neumann boundary conditions, you modified the source term of the elliptic equation?

Starting a guess value = 0 everywhere is not optimal... you can start this way for the first dt, then is better to start from the field at the previous time step.

However, if the iterative solver does not converge for omega = 1 (simple Gauss-Seidel) then you implemented in a wrong way the FTS.

Thanyou for the reply sir.

STEPS about how I solve for poisson equation

1) I first solve the poisson equation for pressure for the INTERNAL scalar grid points.

2) Then , I apply the zero pressure gradient at the boundaries. Now I have have the pressure at ALL the scalar grid points

3) Then, I update the velocities at INTERNAL u and v control volumes ( I already implement dirichlet B.C ) using

u_new(i,j) = u_star(i,j) + (dt/rho/dx) * ( p(i,j) - p(i+1,j) );

v_new(i,j) = v_star(i,j) + (dt/rho/dy) * ( p(i,j) - p(i,j+1) );

4) Then I calculate mass source

rho*dy*( u_new(i,j) - u_new(i-1,j) ) + rho*dx*( v_new(i,j) - v_new(i,j-1) )

at all the INTERNAL scalar C.V

5) then I say,

u_star = u_new;
v_star = v_new;

I do steps ( 1 - 5 ) till the mass source is less than 1e-7.

so , yes in the step 4 I do modify the source term.

And yes I initialize pressure as zero only during the first dt, after that I use the pressures at the old time step.

My code converges when w_opt = 1 , but the problem is that it takes lot of time ( as I mentioned in my previous reply) to reach to the steady state. So, thats what the problem is and so I decided to use SOR but Its not providing any help.

[FMDenaro]

but converging towards the steady state is not an issue due to SOR ... it is involved the physics of the flow and the time for a steady state (until it existes) is quite proportional to the Re number.

furthermore, you cannot solve step 1) without simoultaneously fixing BC.s 2)... Thus I wonder what BC.s are you setting in step 1)

[mihirmakwana6]

Originally Posted by FMDenaro

but converging towards the steady state is not an issue due to SOR ... it is involved the physics of the flow and the time for a steady state (until it existes) is quite proportional to the Re number.

furthermore, you cannot solve step 1) without simoultaneously fixing BC.s 2)... Thus I wonder what BC.s are you setting in step 1)

ok.
sir, I am not setting any B.C when solving step1 , the reason being that I am solving the poisson equation at INTERNAL points and then setting the values at the boundaries such that zero pressure gradient is obtained.


This is the implementation of B.C

p(1,: ) = p(2,: ); % left wall

p(imax,: ) = p(imax-1,: ); % right wall

p( :,1 ) = p( :,2); % bottom wall

p( :,jmax ) = p( :,jmax-1); % top wall

where imax and jmax are maximum points in x and y direction.


Thus after step 2 , I have pressure values at all grid points.

Now, say that after 1st ITERATION , i reach step 1, then , of the entire imax X jmax matrix of pressure , only the matrix of size

i = 2:imax-1 and j = 2:jmax-1 is changed in step 1 and

then in step 2 , the remaining elements change because of the zero gradient B.C.

this happen until the mass source is less than 1e-7.

Is there a problem with this ?

[FMDenaro]

Originally Posted by mihirmakwana6

ok.
sir, I am not setting any B.C when solving step1 , the reason being that I am solving the poisson equation at INTERNAL points and then setting the values at the boundaries such that zero pressure gradient is obtained.


This is the implementation of B.C

p(1,: ) = p(2,: ); % left wall

p(imax,: ) = p(imax-1,: ); % right wall

p( :,1 ) = p( :,2); % bottom wall

p( :,jmax ) = p( :,jmax-1); % top wall

where imax and jmax are maximum points in x and y direction.


Thus after step 2 , I have pressure values at all grid points.

Now, say that after 1st ITERATION , i reach step 1, then , of the entire imax X jmax matrix of pressure , only the matrix of size

i = 2:imax-1 and j = 2:jmax-1 is changed in step 1 and

then in step 2 , the remaining elements change because of the zero gradient B.C.

this happen until the mass source is less than 1e-7.

Is there a problem with this ?

you cannot set zero pressure gradient AFTER the solution of the elliptic equation, it is part of the problem!
By a numerical point of view, you have in the matrix some values therefore during the SOR iteration for the inner points you are implicitly running a fictious Dirichlet problem.

I suggest you studyng better how to solve Neumann problem. Furthermore, after each time step check the continuity equation in each computational cell.

[mprinkey]

By not including the zero-gradient boundary condition (and the necessary single fixed point somewhere in the domain) you are NOT solving the Poisson system. You are solving a portion of the Poisson equation (over the internal cells/points) and then explicitly updating the boundary values. So, no matter if you use a super-accurate/efficient multigrid or FFT solver for the internal part, the FULL SYSTEM is only coupled by some block Gauss Siedel...one block being the internal cells and the other block being the boundary cells.

As Filippo has noted, this is not the way to do this. You need to put all of the pressure unknowns inside the matrix that goes to the Poisson solver. Internal cells get the normal Laplacian entries. But boundary conditions need to be included implicitly. For zero-gradient conditions and low-order schemes, it is easy: P_wall = P_oneCellFromWall. You can keep the same number of unknowns that you currently have, but you have to modify the Laplacian entries for the cells next to the wall. Here is a simple example in 1D.

At the wall for the zero-derivative condition:

P[0] = P[1].

The normal Laplacian entry for cell 1 is something like this (with R[] as the right-hand side source term):

P[0] - 2.0*P[1] + P[2] = R[1]

To include the BC implicitly, you need to modify that cell 1 entry using the first equation:

P[1] - 2.0*P[1] + P[2] = R[1] ... or just

-P[1] + P[2] = R[1]

Of course, other matrix entries away from the boundary remain unchanged:

P[1] - 2.0*P[2] + P[3] = R[2]

When you make that change and implicitly apply the zero-gradient condition on all boundaries, the matrix for pressure becomes indeterminate. So, you need to pick one entry in the matrix and replace it with just this form. You can use a fixed value of P = 0

P[fixedPoint] = P_fixedValue

That will make the system determinate and any correctly coded sparse linear equation solver should solve it robustly and efficiently.

Then, if you like, you can go back and use the boundary condition formula to update the wall values:

P[0] = P[1]

The key difference between this and what you seem to be doing is that modification of the matrix entries for cells touching the boundary. Without that implicit treatment, your code will never converge efficiently and you will be lucky if it will converge at all.

[FMDenaro]

Originally Posted by mprinkey

By not including the zero-gradient boundary condition (and the necessary single fixed point somewhere in the domain) you are NOT solving the Poisson system. You are solving a portion of the Poisson equation (over the internal cells/points) and then explicitly updating the boundary values. So, no matter if you use a super-accurate/efficient multigrid or FFT solver for the internal part, the FULL SYSTEM is only coupled by some block Gauss Siedel...one block being the internal cells and the other block being the boundary cells.

As Filippo has noted, this is not the way to do this. You need to put all of the pressure unknowns inside the matrix that goes to the Poisson solver. Internal cells get the normal Laplacian entries. But boundary conditions need to be included implicitly. For zero-gradient conditions and low-order schemes, it is easy: P_wall = P_oneCellFromWall. You can keep the same number of unknowns that you currently have, but you have to modify the Laplacian entries for the cells next to the wall. Here is a simple example in 1D.

At the wall for the zero-derivative condition:

P[0] = P[1].

The normal Laplacian entry for cell 1 is something like this (with R[] as the right-hand side source term):

P[0] - 2.0*P[1] + P[2] = R[1]

To include the BC implicitly, you need to modify that cell 1 entry using the first equation:

P[1] - 2.0*P[1] + P[2] = R[1] ... or just

-P[1] + P[2] = R[1]

Of course, other matrix entries away from the boundary remain unchanged:

P[1] - 2.0*P[2] + P[3] = R[2]

When you make that change and implicitly apply the zero-gradient condition on all boundaries, the matrix for pressure becomes indeterminate. So, you need to pick one entry in the matrix and replace it with just this form. You can use a fixed value of P = 0

P[fixedPoint] = P_fixedValue

That will make the system determinate and any correctly coded sparse linear equation solver should solve it robustly and efficiently.

Then, if you like, you can go back and use the boundary condition formula to update the wall values:

P[0] = P[1]

The key difference between this and what you seem to be doing is that modification of the matrix entries for cells touching the boundary. Without that implicit treatment, your code will never converge efficiently and you will be lucky if it will converge at all.

Good explanation

I would observe:
1) the RHS of the pressure equation must be properly modified in such a way that the compatibility relation is satisfied. That means to have the volume integral of this term vanishing.

2) fixing a value is not strictly necessary. The solver implicitly fix some value during iterations and you have convergence. It is my experience that fixing a-priori a value will slow-down the convergence.

[mprinkey]

Originally Posted by FMDenaro

Good explanation

I would observe:
1) the RHS of the pressure equation must be properly modified in such a way that the compatibility relation is satisfied. That means to have the volume integral of this term vanishing.

2) fixing a value is not strictly necessary. The solver implicitly fix some value during iterations and you have convergence. It is my experience that fixing a-priori a value will slow-down the convergence.

(1) Right. I should not assume that the div(U) is being done correctly, but we walk...then we run, right?

(2) I don't disagree with your point, but that seems to be a case of the linear solver being smarter than the user. dP/dn =0 and Laplacian(P) = R does lead to an indeterminate system. Well-designed linear equation solvers may often have a better way to impose determinacy than fixing a single point, but a student-written SOR or CG solver likely won't.

[mihirmakwana6]

Originally Posted by mprinkey

By not including the zero-gradient boundary condition (and the necessary single fixed point somewhere in the domain) you are NOT solving the Poisson system. You are solving a portion of the Poisson equation (over the internal cells/points) and then explicitly updating the boundary values. So, no matter if you use a super-accurate/efficient multigrid or FFT solver for the internal part, the FULL SYSTEM is only coupled by some block Gauss Siedel...one block being the internal cells and the other block being the boundary cells.

As Filippo has noted, this is not the way to do this. You need to put all of the pressure unknowns inside the matrix that goes to the Poisson solver. Internal cells get the normal Laplacian entries. But boundary conditions need to be included implicitly. For zero-gradient conditions and low-order schemes, it is easy: P_wall = P_oneCellFromWall. You can keep the same number of unknowns that you currently have, but you have to modify the Laplacian entries for the cells next to the wall. Here is a simple example in 1D.

At the wall for the zero-derivative condition:

P[0] = P[1].

The normal Laplacian entry for cell 1 is something like this (with R[] as the right-hand side source term):

P[0] - 2.0*P[1] + P[2] = R[1]

To include the BC implicitly, you need to modify that cell 1 entry using the first equation:

P[1] - 2.0*P[1] + P[2] = R[1] ... or just

-P[1] + P[2] = R[1]

Of course, other matrix entries away from the boundary remain unchanged:

P[1] - 2.0*P[2] + P[3] = R[2]

When you make that change and implicitly apply the zero-gradient condition on all boundaries, the matrix for pressure becomes indeterminate. So, you need to pick one entry in the matrix and replace it with just this form. You can use a fixed value of P = 0

P[fixedPoint] = P_fixedValue

That will make the system determinate and any correctly coded sparse linear equation solver should solve it robustly and efficiently.

Then, if you like, you can go back and use the boundary condition formula to update the wall values:

P[0] = P[1]

The key difference between this and what you seem to be doing is that modification of the matrix entries for cells touching the boundary. Without that implicit treatment, your code will never converge efficiently and you will be lucky if it will converge at all.

thankyou Michael Prinkey for the reply.


1)

I have modified the poisson equation by incorporating the B.C implicity. thus for the 2D domain , I have mofidfied all the equations for the points near the boundary

i.e j = 2 AND j = jmax - 1 , i = 2 TO imax-1 and

i=2 AND imax - 1 , j = 3 TO jmax - 2

where imax , jmax are maximum points in x and y directions



But, I have not understood WHY is there a need to use

P[fixedPoint] = P_fixedValue

Also, which point in the domain should I select as a fixedpoint.
Also, at what point in the code should I incorporate the above condition ?

Is there a link/paper that discusses this ?

2)

My loop for pressure runs till mass source is GREATER than 1e-7. so , when the mass source becomes LESS than 1e-7 , it itself justifies that the incompressibility constraint is satisfied at each of the two RK2 steps.

[mihirmakwana6]

Originally Posted by FMDenaro

Good explanation

I would observe:
1) the RHS of the pressure equation must be properly modified in such a way that the compatibility relation is satisfied. That means to have the volume integral of this term vanishing.

2) fixing a value is not strictly necessary. The solver implicitly fix some value during iterations and you have convergence. It is my experience that fixing a-priori a value will slow-down the convergence.

sir, by compatibility relation do you mean the incompressibilty constraint ,

div(rho*u) = 0

[FMDenaro]

Originally Posted by mihirmakwana6

thankyou Michael Prinkey for the reply.


1)

I have modified the poisson equation by incorporating the B.C implicity. thus for the 2D domain , I have mofidfied all the equations for the points near the boundary

i.e j = 2 AND j = jmax - 1 , i = 2 TO imax-1 and

i=2 AND imax - 1 , j = 3 TO jmax - 2

where imax , jmax are maximum points in x and y directions



But, I have not understood WHY is there a need to use

P[fixedPoint] = P_fixedValue

Also, which point in the domain should I select as a fixedpoint.
Also, at what point in the code should I incorporate the above condition ?

Is there a link/paper that discusses this ?

2)

My loop for pressure runs till mass source is GREATER than 1e-7. so , when the mass source becomes LESS than 1e-7 , it itself justifies that the incompressibility constraint is satisfied at each of the two RK2 steps.

1) no matter where you fix the value, however, I suggest first trying a boundary point

2) No, convergence of the pressure equation at certain residual level does not imply you have satisfied continuity at the same level. What is more, you can have tha case of some convergence on the pressure equation that can be wrong and does not satisfy the continuity equation in each cell. You must definitely check the continuity constraint after that the velocity (v* - dt Grad p ) is computed.

[FMDenaro]

Originally Posted by mihirmakwana6

sir, by compatibility relation do you mean the incompressibilty constraint ,

div(rho*u) = 0

No ....

From the pressure equation you have

Div Grad p = Div q

and if homogeneous Neumann BC.s are used, the equation above implies:

Int [S] dp/dn dS = Int [S] n.q dS = 0

that is the compatibility equation

[mprinkey]

Originally Posted by mihirmakwana6

But, I have not understood WHY is there a need to use

P[fixedPoint] = P_fixedValue

There are mathematical ways to look at it, but it is really a manifestation of the physics. Pressure in incompressible flow is a gauge variable. You can see that because the incompressible momentum equation only has a grad(P) term. So, you could replace P with P + C, where C is a constant value and you will get exactly the same answer...grad(P + C) = grad(P).

The Poisson equation for P (with zero-gradient BCs) is thus indeterminate. Because any solution P[i] could be replaced with P[i] + C. Iterative solutions will not converge because the "constant" value can drift up and down. There is no right value for C.

Now, if you have a fixed pressure boundary condition, this issue goes away because one or more pressure values are fixed and that forces the C to a single value. In the case of all zero-gradient conditions, you have to explicitly fix one value to lock down the C and prevent solution drift.

[mihirmakwana6]

Originally Posted by mprinkey

By not including the zero-gradient boundary condition (and the necessary single fixed point somewhere in the domain) you are NOT solving the Poisson system. You are solving a portion of the Poisson equation (over the internal cells/points) and then explicitly updating the boundary values. So, no matter if you use a super-accurate/efficient multigrid or FFT solver for the internal part, the FULL SYSTEM is only coupled by some block Gauss Siedel...one block being the internal cells and the other block being the boundary cells.

As Filippo has noted, this is not the way to do this. You need to put all of the pressure unknowns inside the matrix that goes to the Poisson solver. Internal cells get the normal Laplacian entries. But boundary conditions need to be included implicitly. For zero-gradient conditions and low-order schemes, it is easy: P_wall = P_oneCellFromWall. You can keep the same number of unknowns that you currently have, but you have to modify the Laplacian entries for the cells next to the wall. Here is a simple example in 1D.

At the wall for the zero-derivative condition:

P[0] = P[1].

The normal Laplacian entry for cell 1 is something like this (with R[] as the right-hand side source term):

P[0] - 2.0*P[1] + P[2] = R[1]

To include the BC implicitly, you need to modify that cell 1 entry using the first equation:

P[1] - 2.0*P[1] + P[2] = R[1] ... or just

-P[1] + P[2] = R[1]

Of course, other matrix entries away from the boundary remain unchanged:

P[1] - 2.0*P[2] + P[3] = R[2]

When you make that change and implicitly apply the zero-gradient condition on all boundaries, the matrix for pressure becomes indeterminate. So, you need to pick one entry in the matrix and replace it with just this form. You can use a fixed value of P = 0

P[fixedPoint] = P_fixedValue

That will make the system determinate and any correctly coded sparse linear equation solver should solve it robustly and efficiently.

Then, if you like, you can go back and use the boundary condition formula to update the wall values:

P[0] = P[1]

The key difference between this and what you seem to be doing is that modification of the matrix entries for cells touching the boundary. Without that implicit treatment, your code will never converge efficiently and you will be lucky if it will converge at all.

1)

sir, just to clarify that what I am doing is right or not , regarding implementing the condition

P[fixedPoint] = P_fixedValue

say , in the 1D example you mentioned i use P[0] = 0 i.e i take the wall point i = 0 as the fixed point .

Now, normal equation is

P[0] - 2.0*P[1] + P[2] = R[1]

and Now since

P[0] = 0

hence the modified equation becomes

- 2.0*P[1] + P[2] = R[1]

so I have to solve the above equation and NOT the equation that I obtain for point i = 1 by applying Neumann B.C i.e

-P[1] + P[2] = R[1]

Is this right?


2) I am using Red-Black Successive Over-Relaxation (SOR) , where red represents the even grid points and black denotes the odd grid points. I am initially solving at only red point and then the black one. So, is it fine or will have any adverse effect on the solution.

[mihirmakwana6]

Originally Posted by FMDenaro

1) no matter where you fix the value, however, I suggest first trying a boundary point

2) No, convergence of the pressure equation at certain residual level does not imply you have satisfied continuity at the same level. What is more, you can have tha case of some convergence on the pressure equation that can be wrong and does not satisfy the continuity equation in each cell. You must definitely check the continuity constraint after that the velocity (v* - dt Grad p ) is computed.

sir, I did check the mass source term

rho*dy*( u_new(i,j) - u_new(i-1,j) ) + rho*dx*( v_new(i,j) - v_new(i,j-1) )

for each scalar C.V and the maximum value is less than 1e-7 .

which justifies that the pressure is tuned in each scalar cell such that mass continuity is satisfied

[mprinkey]

Originally Posted by mihirmakwana6

1)

sir, just to clarify that what I am doing is right or not , regarding implementing the condition

P[fixedPoint] = P_fixedValue

say , in the 1D example you mentioned i use P[0] = 0 i.e i take the wall point i = 0 as the fixed point .

Now, normal equation is

P[0] - 2.0*P[1] + P[2] = R[1]

and Now since

P[0] = 0

hence the modified equation becomes

- 2.0*P[1] + P[2] = R[1]
...

2) I am using Red-Black Successive Over-Relaxation (SOR) , where red represents the even grid points and black denotes the odd grid points. I am initially solving at only red point and then the black one. So, is it fine or will have any adverse effect on the solution.

(1) That is correct. At that fixed point, you are not applying the zero-gradient condition. You are just setting it to 0.

(2) The linear equation solver should not be an issue. When the matrix is right, it is right. The process you use to solve it is immaterial--assuming that the linear equation solver is coded correctly.

[mihirmakwana6]

I found the bug in my code. Actually I was modifying the source term in the poisson solver which actually needs to be fixed till the pressure loop runs. Now, My code runs well and I get my results quite fast.

Also, if I use

P[fixedPoint] = P_fixedValue

than it takes more time then the case when i don't use it.

Thankyou Filippo Maria Denaro sir and Michael Prinkey sir for all your replies

- Mihir