NaN causing assertion failure in PsyRhovFnTdbWPb_fast

[JasonGlazer]

The user file stops at an assertion failure when running in develop because a NaN (shown in Windows as a -1.#IND000000000000) was passed to PsyRhovFnTdbWPb_fast. The file was provided by a user as an 8.1 file that crashed. Tested on Windows.

• [X] Defect file added see: EnergyPlusDevSupport/tree/master/DefectFiles/5000s/5520
• [x] Ticket added to Pivotal for defect (development team task)
• [X] Pull request created (the pull request will have additional tasks related to reviewing changes that fix this defect)

[rraustad]

It would be nice to be able so set a compiler switch to break on this error.

[mjwitte]

@JasonGlazer posted this recently https://github.com/NREL/EnergyPlus/wiki/Debugging-Tips

[rraustad]

line 3922 of ConvectionCoefficients: HcIn( SurfNum ) = std::pow( std::pow( HcIn( SurfNum ), 3.2 ) + std::pow( Hf, 3.2 ), 1.0 / 3.2 ); SurfNum = 5 and HcIn is allocated to 12 HcIn(5) = 0.72122 Hf = -7.337

line 331 of CrossVentMgr: CalcDetailedHcInForDVModel( SurfNum, TempSurfIn, CVHcIn, Urec );

CalcDetailedHcInForDVModel(
int const SurfNum, // surface number for which coefficients are being calculated
Array1S< Real64 > const SurfaceTemperatures, // Temperature of surfaces for eval of HcIn
Array1S< Real64 > HcIn, // Interior Convection Coeff Array
Optional< Array1S< Real64 > const > Vhc // Velocity array for forced convection coeff calculation

Hf = 4.3 * Vhc()( Surface( SurfNum ).Zone ); Vhc is -1.706 which was passed in from Urec in the call above.

I think this is coming from cross vent model with something to do with the jet flow. I can't decipher.

[JasonGlazer]

@rraustad I believe that the line

https://github.com/NREL/EnergyPlus/blob/develop/src/EnergyPlus/ConvectionCoefficients.cc#L3906

Hf = 4.3 * Vhc()( Surface( SurfNum ).Zone );

Is converting from a velocity to a convection coefficient but since Vhc is coming from the room air models it seems like it can be both positive and negative values (I guess am not familiar with that stuff). The convection coefficient should only be dependent on the velocity independent of direction so I think this would be valid:

Hf = 4.3 * abs(Vhc()( Surface( SurfNum ).Zone ));

This takes care of the NaN. Does that seem reasonable?

[mjwitte]

Vhc is an array? How does that work here?

[rraustad]

I agree that line is where this NaN issue started, although I don't know if a negative number for Vhc is valid. I'd first suggest not wrapping the Surface.Zone inside the abs(). I'd then like to know if the example file(s) changes. I'd then like to know when/why Vhc is negative.

Maybe @EnergyArchmage could help with answers?

[JasonGlazer]

Optional< Array1S< Real64 > const > Vhc // Velocity array for forced convection coeff calculation.

The Vhc()(index) syntax seem weird by it seems to just be accessing that item in the array.

[JasonGlazer]

@rraustad and @EnergyArchmage It looks like the first negative value is being calculated at:

https://github.com/NREL/EnergyPlus/blob/develop/src/EnergyPlus/CrossVentMgr.cc#L763

CVJetRecFlows( Ctd, ZoneNum ).Vjet = CVJetRecFlows( Ctd, ZoneNum ).Uin * std::sqrt( CVJetRecFlows( Ctd, ZoneNum ).Area ) * 6.3 * std::log( Dstar( ZoneNum ) / ( 6.0 * std::sqrt( CVJetRecFlows( Ctd, ZoneNum ).Area ) ) ) / Dstar( ZoneNum );

The problem is that this part:

Dstar( ZoneNum ) / ( 6.0 * std::sqrt( CVJetRecFlows( Ctd, ZoneNum ).Area ) )

a small value and taking the std::log of it makes it negative. This line of code is directly discussed in the Engineering Reference as equation 147 on page 492 (EnergyPlus 8.4 release). My guess is that the correlation just doesn't make too much sense for the values used. In this case Dstar = 14.04 and the Area is 31.1. Unadjusted Droom is 13.63.

Suggestions?

[rraustad]

Vjet shouldn't be negative. Not sure what to do here when that happens.

[EnergyArchmage]

Lets see.

I liked the original suggestion to do abs(Vhc).

I doesn't necessarily suprise me that a room air flow velecity could be negative to indicate a backflow directions, but I'd have to dig into this model (from G.) to be sure if that is okay,

That does seem a likely candidate for an numerical problem. I would do an if trap for " Dstar( ZoneNum ) / ( 6.0 * std::sqrt( CVJetRecFlows( Ctd, ZoneNum ).Area ) )" >= 1.0 before that log expression. Then the else would just replace the enitre log term with 1.0.

[rraustad]

So then 6*sqrt(Area)=Dstar is a limit. Geometrically this seems like a limit on the width of a space given the depth. I did ask myself which space this was that had a depth = 14m while area = 31m, thought it might be a hallway with a 2m width.

[JasonGlazer]

Thank you @rraustad and @EnergyArchmage I have added a check to the routine that does this and it no longer crashes.