[FLASH-BUGS] cylindrical coordinates, avisco

[Tomasz Plewa]

Matt and All -

It's fixed in the next version. Use the following:

!-------  (r_cyl, z)
!
     else if (igeomx == 1  .and.  igeomy == 0) then 
!
        if (xyzswp == 1) then 

           dytb = 0.5e0 / (ytop - ybot)

           do i = nzni, nznf
              dxtb = x(i)**2 - x(i-1)**2
              if ( dxtb /= 0.e0 ) then
                 avis(i) = (x(i) * u(i) - x(i-1) * u(i-1)) * 2.e0 / dxtb   &
                          +(uttp(i) + uttp(i-1) - utbt(i) - utbt(i-1))     &
                          *dytb
                 avis(i) = - cvisc * avis(i) * (x(i) - x(i-1))
              end if
           enddo

        else

           do i = nzni, nznf
              dxtb = xtop**2 - xbot**2
              if ( dxtb /= 0.e0 ) then
                 avis(i) = (u(i) - u(i-1)) / (x(i) - x(i-1))               &
                          +( xtop * (uttp(i) + uttp(i-1))                  &
                            -xbot * (utbt(i) + utbt(i-1))) / dxtb
                 avis(i) = - cvisc * avis(i) * (x(i) - x(i-1))
              end if
           enddo

        end if

Tomek
---
On Fri, Sep 26, 2003 at 03:47:23PM -0400, Sean Matt wrote:
> 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

--