[FLASH-USERS] Typo in unsplit_rad Hydro_computeDT

[Klaus Weide]

On Mon, 8 Feb 2016, Mark Richardson wrote:

> >   I'm pretty sure there's a bug in the unsplit_Rad Hydro_computeDT:
> > source/physics/Hydro/HydroMain/unsplit_rad/Hydro_computeDt.F90
> >
> > Line 241 reads:
> > 241:                 if (NDIM > 2) dt_ltemp =
> > max(dt_ltemp,(abs(U(VELZ_VAR,i,j,k)-uzgrid(j))+sqrt(cfz2))*delzinv)
> >
> > But should be:
> > 241:                 if (NDIM > 2) dt_ltemp =
> > max(dt_ltemp,(abs(U(VELZ_VAR,i,j,k)-uzgrid(k))+sqrt(cfz2))*delzinv)

Hello Mark,

Agreed on the bug. It also affects the version of Hydro_computeDt in 
source/physics/Hydro/HydroMain/unsplit/Hydro_computeDt.F90.
This will be fixed in the next release.

As long as as the arrays uxgrid/uygrid/uzgrid are zero-filled, which is 
always the case in the usual calls Hydro_computeDt by other FLASH code,
this bug had the potential to cause problems only when the following 
conditions were satisfied:
 (1) 3D setup with unsplit Hydro or MHD;
 (2) non-cubic (logical) block sizes - specifically, NZB > NYB.

Thank you for reporting the bug.

> > However, I was getting the error:
> > ***********************************************************
> >   Warning: The initial timestep is too large.
> >     initial timestep =   1.000000000000000E-010
> >     CFL timestep     =   5.215016855876254E-303
> >   Resetting dtinit to dr_tstepSlowStartFactor*dtcfl.
> >  ***********************************************************
> > both with and without the fix above. With the fix, I dug deeper into why
> > the dt was going to zero, and found that it related to guard cells. I then
> > also initialized my guard cells in Simulation_initBlock and the small dt
> > problem went away.

This is a bit mysterious. Hydro_computeDt should not depend on any values 
in the guard cells of the current block at all, unless BOTH of the 
following are true: (A) A local CFL factor is being used (i.e., CFL_VAR
is defined); and (B) you have changed the Fortran parameter hy_cflStencil
in Hydro_data from 0 to something > 0. Have you made any modifications?

Where and how did you fill guard cells at initialization?

The following assumes that you are using PARAMESH as your Grid 
implementation (not UG, not PARAMESH2).

In general, initializing guard cells in your Simulation_initBlock should 
not be necessary (but cannot hurt either). The logic for calling 
Simulation_initBlock and then calling Eos and potentially filling guard 
cells (indirectly, via Grid_markRefineDerefine) is in gr_expandDomain.F90,
q.v.  The first calls to Hydro_computeDt come from Driver_verifyInitDt.
It is indeed possible that at the point when Driver_verifyInitDt calls
Hydro_computeDt, guard cells have NOT been properly filled or updated,
although often they will have been updated as an effect of what happened
in gr_expandDomain. But again, Hydro_computedDt should not normally need
them (except when (A)&&(B) apply, see above); that's why the 
    call Grid_fillGuardCells
in  Driver_verifyInitDt.F90 is commented out.

> > Oddly, a collaborator took my SimulationMain/Code/ and ran it, and got the
> > small dt error again, even with the guardCells initialized. I forgot to
> > give him the above fix to Hydro_computeDT. He made that fix, and now his
> > small dt problem is gone as well.
> >
> > I'm using the hypre unsplit solver, with radiation and mgd on. Any
> > thoughts on why an intiial FillGuardCells call isn't fixing the problem?