[FLASH-BUGS] cylindrical coordinates, avisco

[Sean Matt]

Thanks, Tomek, for the quick response.  Could you please tell me the 
proper way to declare/define nzni, and nznf (which are not used in the 
original 2.3 release version of avisco.F90)?  Thanks.

	-Sean

On Fri, 26 Sep 2003, Tomasz Plewa wrote:

> 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
> 
> 

-- 

_______________________________________________________________________
Sean P. Matt				     CITA National Fellow
phone: (905) 525-9140 x 23189		     Physics & Astronomy Dept.
http://www.physics.mcmaster.ca/~matt	     McMaster University
-----------------------------------------------------------------------