Bug found in the coordinate class for schwarzschild geometry
Closed this issue · 1 comments
Dear developers,
upon using Athena I came across a bug that I would like to turn you attention to.
On line 356 and 357 of src/coordinates/schwarzschild.cpp
athena/src/coordinates/schwarzschild.cpp
Line 357 in 185473d
one can read
// \Delta W = r \Delta\theta
dx2(i) = x1v(i) * dx1f(j);
when one should actually reads
// \Delta W = r \Delta\theta
dx2(i) = x1v(i) * dx2f(j);
as the previous one is actually equivalent to dx2 = r dr, which is numerically incorrect and also led me to weird problems in running the code, specially when my meshblocks were not defined as a cube. This came to be as dx1f has length of max(i) and dx2f has length of max(j). In the cube meshblocks as max(i)=max(j) this was no problem, but for generic meshblock structures I was systematically get dt=0, as dt ~ min(dx1, dx2,dx3) and dx2 was set to zero as for j >max(i) dx1f(j) was not not defined.
Agreed, this is indeed a bug. Thanks for bringing it to our attention!