Error in AddFeedbackSphere

[Elizabeth Tasker]

Hi,

I just wanted to check this is really an error before I correct.

In Grid_AddFeedbackSphere, line 821 under "STROEMGREN SPHERE FROM A
POP III STAR (WHALEN ET AL. 2004)" there is:

if (GENum >= 0)
BaryonField[GENum][index] = EjectaThermalEnergy*ramp;

But it should be:

if (GENum >= 0 && DualEnergyFormalism)
BaryonField[GENum][index] = EjectaThermalEnergy*ramp;

Since otherwise, if DualEnergyFormalism = 0, then GENum = 0 which is
the same number as DensNum.

In an unrelated test problem where I made this mistake, I got
unexpected supernovae-esque explosions. I rather regret tracing it to
a bug and not just declaring they were violent cloud collisions....

(Grid_PutSinkRestartInitialize.C also has this error, but only in a
commented-out section).

Elizabeth

[Matthew Turk]

Hi Elizabeth,

On Wed, Feb 22, 2012 at 10:25 AM, Elizabeth Tasker
<tas...@astro1.sci.hokudai.ac.jp> wrote:
> Hi,
>
> I just wanted to check this is really an error before I correct.
>
> In Grid_AddFeedbackSphere, line 821 under "STROEMGREN SPHERE FROM A
> POP III STAR (WHALEN ET AL. 2004)" there is:
>
>  if (GENum >= 0)
>              BaryonField[GENum][index] = EjectaThermalEnergy*ramp;
>
> But it should be:
>
>  if (GENum >= 0 && DualEnergyFormalism)
>              BaryonField[GENum][index] = EjectaThermalEnergy*ramp;

I agree. In Grid_IdentifyPhysicalQuantities.C, the variables are
initialized to 0:

DensNum = GENum = Vel1Num = Vel2Num = Vel3Num = TENum = 0;

Following this, GENum is set to FindField, which *can* be less than
zero, but *only* if DualEnergyFormalism == TRUE. (A change I would
also suggest to your proposed change above).

/* Find gas energy, if possible. */

if (DualEnergyFormalism == TRUE)
if ((GENum = FindField(InternalEnergy, FieldType,
NumberOfBaryonFields)) < 0) {
ENZO_FAIL("Cannot find gas energy.");
}

This probably did not ever come up because to my recollection
STROEMGREN is only turned on for DEF simulations. Very good catch.

I would recommend submitting a pull request with this. For the
unrelated test problem, is it one that would be valuable to add to the
test suite, to catch regressions like this?

-Matt

>
> Since otherwise, if DualEnergyFormalism = 0, then GENum = 0 which is
> the same number as DensNum.
>
> In an unrelated test problem where I made this mistake, I got
> unexpected supernovae-esque explosions. I rather regret tracing it to
> a bug and not just declaring they were violent cloud collisions....
>
> (Grid_PutSinkRestartInitialize.C also has this error, but only in a
> commented-out section).
>
> Elizabeth
>

> --
> You received this message because you are subscribed to the Google Groups "enzo-dev" group.
> To post to this group, send email to enzo...@googlegroups.com.
> To unsubscribe from this group, send email to enzo-dev+u...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/enzo-dev?hl=en.
>

[Elizabeth Tasker]

Hi Matt,

>> But it should be:
>>
>>  if (GENum >= 0 && DualEnergyFormalism)
>>              BaryonField[GENum][index] = EjectaThermalEnergy*ramp;
>
> I agree.  In Grid_IdentifyPhysicalQuantities.C, the variables are
> initialized to 0:
>
>  DensNum = GENum = Vel1Num = Vel2Num = Vel3Num = TENum = 0;
>
> Following this, GENum is set to FindField, which *can* be less than
> zero, but *only* if DualEnergyFormalism == TRUE.  (A change I would
> also suggest to your proposed change above).

So maybe change Grid_IdentifyPhysicalQuantities to:

GENum = FindField(InternalEnergy, FieldType,
NumberOfBaryonFields);

if (DualEnergyFormalism == TRUE && GENum < 0)

ENZO_FAIL("Cannot find gas energy.");

Since it's somewhat asking for disaster to have GENum set to 0.

>  For the
> unrelated test problem, is it one that would be valuable to add to the
> test suite, to catch regressions like this?

Actually, I used the O'Shea definition of "Test Problem" to mean "An
Excellent Simulation of Notable Physical Worth that just happens not
to be Cosmological". It was a restart of an isolated galaxy disc
(similar to SupernovaeRestartInitialize.C), so I don't think on it's
own it's that useful.

That said, I do wonder if there is a way of sanity checking the field
numbers since there are now so many, especially since I found an
overlap in the typdef.h the other day. Since they're all set
automatically (or they should be) by routines like
IdentifyPhysicalQuantities, FindField and AddFields, a test problem
might be the way to check those all play nicely together. Maybe just
one that queried all possible fields and ensured the answer back was
sensible. I'm not sure how to make that non-trivial.

If you agree with the change to Grid_IdentifyPhysicalQuantities, I'll
make it and send a pull request.

Elizabeth