[FLASH-BUGS] cylindrical coordinates, avisco

Sean Matt matt at physics.mcmaster.ca
Fri Sep 26 14:47:23 CDT 2003


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




More information about the flash-bugs mailing list