Hi All,
We've recently been checking out the 2D cylindrical coordinates in
FLASH for some stellar wind sims, and we have either found a potential
problem or exposed some of our own ignorance. We'm running on Tru64 (OSF1
V5.1) using a Compaq fortran compiler (V5.5-1877-48BBF).
So for our axisymmetric setup, we set xmin = 0.0 (xmax = 1.0) and
the xl_boundary_type = "reflecting". If we are supposed to use a
different boundary type, then our problem is our fault. When we use the
reflecting boundary, the problem initializes properly (and writes the
first output files) but crashes during the first timestep. It crashes
with and FPE on line 116 of avisco.F90. Below is the code lines around
line 116:
do i = 5, nzn4+1
(line 116) avis(i) = (x(i)*u(i) - x(i-1)*u(i-1)) * 2.0 / &
(x(i)**2 - x(i-1)**2) + &
(uttp(i) + uttp(i-1) - utbt(i) - utbt(i-1)) * &
dytb
avis(i) = - cvisc * avis(i) * (x(i) - x(i-1))
enddo
We've tracked down the FPE as a divide by zero. Specifically,
(x(i)**2 - x(i-1)**2) will equal zero for i = 5. This is because, with
the reflecting boundary condition, x(5) = dx/2 and x(4) = -dx/2. The
difference of the squares of these equal and opposite #'s is zero. We can
"get around" this crash by setting our xmin to something like 0.001, but
we think the problem should probably be fixed in the source code.
This actually raises the issue of whether the quantity (x(i)**2 -
x(i-1)**2) is the right formulation for the calculation of the
artificially viscocity. We're not sure why the actual numerical value of
the x coordinate should be relevant, since it is arbitrarily set by the
user (e.g., 201**2-200**2 does not equal 2**2-1**2). Perhaps what is
needed is something like the square of the difference, (x(i) - x(i-1))**2,
but since this is a different quantity, this should be looked at with
care, and we don't know the correct formulation.
-Sean
This archive was generated by hypermail 2b30 : Fri Sep 26 2003 - 14:47:29 CDT