[FLASH-BUGS] FPE in 3D mhd setup

[Timur Linde]

Hi Sean,

We will examine this problem immediately.

Timur

----------------------------------------------------------
Timur Linde                       ASCI Flash Center
Research Scientist                Astronomy & Astrophysics
773-834-3226 (tel)                University of Chicago
773-834-3230 (fax)                5640 South Ellis Ave.
t-linde at uchicago.edu              Chicago, IL 60637
---------------------------------------------------------- 


On Fri, 21 Feb 2003, Sean Matt wrote:

> Hello,
> 
> 	We have been having a few problems getting a new 3D mhd test
> problem setup to run.  We am getting a floating point exception (FPE)
> crash, as well as bad results (before the crash).  We are running on Tru64
> (OSF1 V5.1) using a Compaq fortran compiler (V5.5-1877-48BBF).
> 
> 	The test is the spin down of a magnetic rotating disk via the
> transport of torsional alfven waves along field lines threading the disk.  
> Our setup files are attached to this email; please let me know if you have
> a problem receiving them.  There is an analytical solution for this
> problem (see references listed in Config file), and we think it will be an
> excellent test of the FLASH mhd.
> 	Earlier, we have run the problem using the /hydro module, and with
> no B fields.  The initial setup is a dense, rotating cylinder in a cold
> ambient medium.  The hydro setup seems to work fine, and as the system
> evolves the disk spins apart, as it should (there is nothing keeping it
> together).
> 	Currently, we have switched to the /hydro/mhd module, but 
> initialize magx, magy, magz to zero.  When we run, we get a FPE.  The 
> error occurs after 13 timesteps, in the current setup.  If we initialize 
> magx, magy, magz to 1.0e-10 or to 1.0, the FPE still occurs, though 
> usually at a slightly different timestep.
> 	Using dbx, we've been able to track down the FPE.  It occurs in 
> the routine /source/hydro/mhd/mhd_fluxes.F90 (around line 76) at the 
> statement:
> 
> 	s = sqrt(Um(idn,i)/Up(idn,i-1))
> 
> We found, that at the time of the crash, one of the values of Um and Up
> was positive, while the other was negative.  This results in the attempted
> calculation of the square root of a negative number, hence the FPE.  Um
> and Up are passed in from mhd_sweep.F90 just after being initialized by
> the routine mhd_interpolate.  We gather that Um and Up are
> "time-interpolated data" at the left and right interfaces of the cell.  
> Since the index "idn" refers to the gas density, we're not sure that these
> should ever be negative (though we are only just learning about the 8-wave
> solver).  We were unable to track the problem further.  This is our first 
> major concern.
> 	Second (this may or may not be related to the above), we have
> looked at output data prior to the FPE crash at 13 timesteps.  The data
> doesn't look very good (esp. when compared to the hydro results of the
> same setup at the same time).  In particular, the density in the cylinder,
> which is initially constant, quickly takes on a "mottled" appearance
> (meaning that the values seem to fluxuate in space about the correct
> value).  The pressure in the cylinder also has the wrong value and is also
> mottled.  It is possible that whatever causes the FPE, also causes this
> effect.  We really have no idea.  If this is an actual bug, it would be
> crucial to fix.
> 	We would appreciate any insight you may have about these problems.  
> Please let us know if you require more information.
> 
> 
> 		-Sean
>