Re: [FLASH-BUGS] cylindrical coordinates, avisco

From: Tomasz Plewa (tomek@flash.uchicago.edu)
Date: Fri Sep 26 2003 - 14:51:31 CDT

  • Next message: Tomasz Plewa: "Re: [FLASH-BUGS] cylindrical coordinates, avisco"

    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

    --

    This archive was generated by hypermail 2b30 : Fri Sep 26 2003 - 14:51:33 CDT