[naffrancois]

Hello,

A possible problem may arise from the sub-steps of the runge kutta integrator. Indeed, you apply periodicity on Q, compute the fluxes and update Q1 as Q1=Q+dt*RHS. However this update overrides the ghost cells of Q1 (Q1(0,:,: ) , etc) but you reuse Q1 in the next step for RHS evaluation. Same holds for Q2.

Each conservative vector argument of your function RHS should be applied periodicity as you did it for Q. e.g.:

Q1=Q+dt*RHS
Q1(0:,:, )=Q1(nx,:,: )
... etc
Then RHS(Q1)

Q2=...
Q2(0:,:,: )=Q2(nx,:,: )
... etc
Then RHS(Q2)

If I were you, for debugging purposes, I would switch to a simple Euler time integrator. For 1st order spatial discretization is ok.

The second point to check is the vector RHS. RHS(1), RHS(nx) and RHS(ny) are not defined, but you said indices 1, nx, ny are physical cells.

You should carefully revise your indices, I believe several errors are still there. Ex: Mxp(nx) is not defined because Mxch(nx) is not filled in any value.

Also be careful in your evaluation of dt. It should be inside the temporal loop as your eigenvalues will vary in time.

[selig5576]

Hi,

Thank you very much for that remark. Upon looking at my code more closely I that there are terms that do not have any value at 1, nx, and ny. I have tow thoughts:

1. Extend my arrays so I have (0:nx+2,0:ny+2), so I still have ghost points and do not have to subtract nx-1 so maintain a symmetric array

2. Should I apply boundary conditions at the points of nx, ny, and 1? I'm inclined to do the first option, however would be interested to see your perspective.

Thanks!

[FMDenaro]

Originally Posted by selig5576

Hi,

Thank you very much for that remark. Upon looking at my code more closely I that there are terms that do not have any value at 1, nx, and ny. I have tow thoughts:

1. Extend my arrays so I have (0:nx+2,0:ny+2), so I still have ghost points and do not have to subtract nx-1 so maintain a symmetric array

2. Should I apply boundary conditions at the points of nx, ny, and 1? I'm inclined to do the first option, however would be interested to see your perspective.

Thanks!

Are you using the Fortran compiler with the array bounds check flag?

[selig5576]

Hi Denaro,

I am using the G95 compiler. It will pop up messages if I am trying to access arrays outside the dimension of the array.

[naffrancois]

I do not know compiler options for g95, but I agree, always make sure in development that you use all available checks your compiler offers you, it can save you much much time.

I attached a possible index notation, although you may use whatever you want, the important here is to evaluate quantities where they should be.