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 -----------------------------------------------------------------------
This archive was generated by hypermail 2b30 : Fri Sep 26 2003 - 15:04:04 CDT