[FLASH-USERS] xflash3 and 3D spherical coordinates

[Klaus Weide]

On Fri, 26 Mar 2010, Massimo Gaspari wrote:

> 
> I've tried many tests, changing min max coords, CFL number, NXB vs 
> NYB/NZB. Initial conditions (rho, temp, etc.) seem ok (using visit), but 
> in the first timestep I always receive, with gdb (and -fpe0 enabled), 
> this error in avisco.F90:
>  
> 
> Program received signal SIGFPE, Arithmetic exception.
> avisco (j=5, k=5, avis=..., dirgeom=..., xyzswp=2, nzni=5, nznf=13, nsdim=3, xtop=8.364632500000001e+23, xbot=-2.3767975000000004e+23,
>     ytop=4.4465908125712189e-323, ybot=4.4465908125712189e-323, ylft=4.4465908125712189e-323, yrgt=0, zlft=-0.19634954084936207, zrgt=0.58904862254808621,
>     x=..., xl=..., xzn=..., yzn=..., zzn=..., u=..., uttp=..., utbt=..., utrt=..., utlt=..., cvisc=0.10000000000000001) at avisco.F90:417
> 417                   avis(i) = (sin(x(i)) * u(i) - sin(x(i-1)) * u(i-1)) /        &
> 
> 
> I've print:
> 
> (gdb) print dytb
> $2 = 7.9050503334599447e-323
> (gdb) print dxtb
> $3 = 2.5055296642745135e-72
> (gdb) print dztb
> $4 = 6.9533558064452965e-310

Massimo,

What was your coordinate range for this run; in particular what was the value of xmin?

Could you please also print the values of
dzrl, sinth, i, xzn(i), yzn(i), and zzn(i) when this occurs?

Also, I recommend you make a change in avisco.F90, to split up the long 
expression into several statements, so it should be easier to see which 
subexpression causes the exception. Like so (beginning at line 417):

Replace

              avis(i) = (sin(x(i)) * u(i) - sin(x(i-1)) * u(i-1)) /        &
                   (xzn(j) * sinth * (x(i) - x(i-1)) ) +                   &
                   (xtop**2 * (uttp(i) + uttp(i-1)) -                      &
                   xbot**2 * (utbt(i) + utbt(i-1))  ) * dxtb +             &
                   (utrt(i) + utrt(i-1) - utlt(i) - utlt(i-1)) /           &
                   sinth * dzrl

with

              avis(i) = (sin(x(i)) * u(i) - sin(x(i-1)) * u(i-1)) /        &
                   (xzn(j) * sinth * (x(i) - x(i-1)) )
              avis(i)=avis(i)+ (xtop**2 * (uttp(i) + uttp(i-1)) -            &
                   xbot**2 * (utbt(i) + utbt(i-1))  ) * dxtb +             &
              avis(i)=avis(i)+ (utrt(i) + utrt(i-1) - utlt(i) - utlt(i-1)) / &
                   sinth * dzrl

The error occurs when Hydro is trying to compute and artificial viscosity term for
the Y (==theta) sweep. It seems likely that this is happening when processing
a cell that is very close to x=0.0 (i.e. R=0) and theta=0. It is not clear to
me which subexpression would cause a numerical overflow here, therefore my request
for more data.

---------

As a quick workaround (not: solution), you can try setting runtime parameter cvisc=0.0
(no artificial viscosity). The subroutine avisco should then not even get called.
This should allow you to proceed, maybe until a similar error shows up in another place.



Klaus