Fix floating point exception in FGGasCell when gas contents approach zero

[DheerendraTomar]

When a gas cell is rapidly emptied (e.g. weather balloon burst at high altitude), the gas Contents, Volume, and Pressure can all approach zero. This causes division-by-zero floating point exceptions in three places within FGGasCell::Calculate():

1. GasDensity = GasMass / GasVolume - when GasVolume reaches zero
2. Volume = Contents * R * Temperature / Pressure - when Pressure is zero
3. dVolumeIdeal computation dividing by Pressure and OldPressure

Add guards to handle these zero-denominator cases:

• Return GasDensity = 0 when GasVolume is zero (no gas = no density)
• Set Volume to just BallonetsVolume when Pressure is zero (empty cell)
• Set dVolumeIdeal to 0 when either Pressure is zero (no gas to expand)

Fixes #654

[andgi]

Looks good to me.

Cheers,
Anders Gidenstam

[bcoconni]

I am a bit surprised by this PR:

• Following the commit 464e591, Contents is already guarded from going below $10^{-8}$ so I am struggling to see a scenario where GasVolume == 0.0 ?
  

  jsbsim/src/models/FGGasCell.cpp

  Lines 328 to 329 in aa28a26

  	Contents =
  	max(1E-8, Contents - Pressure * VolumeValved / (R * Temperature));

• Ditto for Pressure which cannot go below the (non-zero) ambient pressure in.Pressure aka AirPressure per the definition of the function max and per the fact that MaxOverpressure is positive.
  

  jsbsim/src/models/FGGasCell.cpp

  Lines 303 to 309 in aa28a26

  	const double IdealPressure =
  	Contents * R * Temperature / (MaxVolume - BallonetsVolume);
  	if (IdealPressure > AirPressure + MaxOverpressure) {
  	Pressure = AirPressure + MaxOverpressure;
  	} else {
  	Pressure = max(IdealPressure, AirPressure);
  	}

Could you provide an example that is actually fixed by your PR ? I mean an example which crashes without your patch and works when your patch is committed.

[AirshipSim]

When a gas cell is rapidly emptied (e.g. weather balloon burst at high altitude), the gas Contents, Volume, and Pressure can all approach zero.

If a cell burst, pressure may approach zero but is at worse Patmosphere, I believe. Also, Patmosphere only reaches close to zero at the Kármán line (100km). A high altitude balloon usually goes up to 40km.

The concept of a gas cell with zero volume is troubling. That would be an empty envelope, a simple piece of fabric. Basically you can put a T-shirt in the air and call it a gas cell with zero volume.

@DheerendraTomar 's case seems indeed exceptional.

[DheerendraTomar]

@bcoconni @AirshipSim Thank you both for the thoughtful feedback.

You're absolutely right that with the current max(1E-8, ...) clamp from commit 464e591, the FPE cannot occur on master. My original PR description was misleading — I should have been clearer about the intent.

The goal of this PR is to replace the temporary fix (464e591, which was explicitly titled "Temporary fix for issue #654") with a proper solution. Here's why the temporary clamp is problematic:

1. It's not physically correct. max(1E-8, ...) means 10⁻⁸ mol of gas can never be vented from the cell, no matter what. The cell can never truly empty. For the gas cell, this is a tiny amount and mostly harmless. But the analogous clamp in FGBallonet is max(1.0, ...) — that's 1 mol of air permanently trapped in the ballonet that can never be removed by any valve operation. That's a noticeable unphysical artifact.

2. The issue is still open. Issue Weather Balloon script crash #654 was left open precisely because 464e591 was a workaround, not a fix. The underlying problem — that the code has no defined behavior for an empty gas cell — was never addressed.

3. Our approach: early return for empty cell. Instead of clamping Contents to an artificial minimum or guarding individual divisions, we allow Contents to reach 0.0 naturally via max(0.0, ...), and then handle the empty state explicitly with an early return. This sets physically correct values:

   • Pressure = AirPressure (an empty envelope equalizes to ambient)
   • Volume = BallonetsVolume (only ballonet volume remains)
   • Mass = 0 (no gas)
   • dVolumeIdeal = 0 (nothing to expand)

   This is analogous to the existing if (Contents > 0) guard already present for the temperature computation at line 288.

@AirshipSim You raise a good point that a zero-volume gas cell is "an empty envelope, a simple piece of fabric." That's exactly right — and that's what happens physically when a balloon bursts. Our early return handles this state cleanly rather than pretending there's always some residual gas.

To reproduce the original crash (before 464e591): revert the max(1E-8, ...) clamp back to max(0.0, ...) in FGGasCell::Calculate() and run scripts/weather-balloon.xml — it will SIGFPE. Our early-return block prevents this crash while also removing the artificial clamp.

[AirshipSim]

weather balloon burst at high altitude

The thing is, what do you actually mean by that? The term 'burst' feels like the cell is utterly destroyed in an instant. In which case, pressure, volume, etc, none of that makes any sense anymore. It sounds like a "crash event" which should be handled out-of-band. For example, the instance of the class representing the gas cell is removed from the "High Altitude Balloon" object. (Note: I have not read the buoyancy code yet, speculative)

Pressure = AirPressure (an empty envelope equalizes to ambient)
Volume = BallonetsVolume (only ballonet volume remains)
Mass = 0 (no gas)

Well, you write Pressure = AirPressure, but isn't Pressure the Pressure of the moles of lifting gas (as in Pressure * V = n(lifting gas) * RT? In which case, Mass(lifting gas) must be different from 0. An empty envelope does not equalize to ambient, it has no volume so its pressure is undefined (P = nRT / 0). This is what I tried to say with my "simple piece of fabric" image. Again, I have not read the code, so I will stop commenting on this PR. I was just trying to understand from the theoretical point of view.

[bcoconni]

Thanks for the clarification, @DheerendraTomar.

However the root of the problem is that we have a mechanism (the valve opening) that leads to the spontaneous creation of vacuum in the envelope. That is the physical meaning of Contents == 0. And this is in violation of the 2nd law of thermodynamics (just as heat moving from cold to hot) so the underlying equations in FGGasCell.cpp are wrong.

The code in FGGasCell.cpp somehow involuntary hides that inconsistency behind the concomitant reduction of the volume as the content vanishes which, when compared with the intuition of a deflating balloon, seems to make sense. But would the envelope be rigid, the equations for the valve opening would also lead to Content == 0.0 meaning that vacuum would be created by the mere fact that the gas density inside the envelope is lower than the density of the ambient gas:

jsbsim/src/models/FGGasCell.cpp

Lines 324 to 327 in 0238e96

	const double DeltaPressure =
	Pressure + CellHeight * g * (AirDensity - GasDensity) - AirPressure;
	const double VolumeValved =
	ValveOpen * ValveCoefficient * DeltaPressure * dt;

That would be the cheapest vacuum pump in the world! Just heat up the air inside an envelope to reduce its density, open a valve and voilà! perfect vacuum has been created in the envelope.

As far as I am concerned your patch does not address the actual problem and is nothing more than a glorified version of my fix, hence your PR remains a temporary fix and as such I am not willing to merge it.

I will not consider the issue #654 fixed until the physics in FGGasCell.cpp are fixed.

[DheerendraTomar]

@bcoconni Thank you — that's an excellent point about the thermodynamics, and I think you've identified the real bug.

I dug into the valve equations and found something interesting. In the manual valving section, GasDensity is computed as:

GasDensity = GasMass / GasVolume = (Contents * M_gas) / (Contents * R * T / P) = M_gas * P / (R * T)

The Contents term cancels out entirely — GasDensity is independent of how much gas remains in the cell. This means DeltaPressure stays positive for any lighter-than-air gas regardless of how much has been vented. The valve model therefore drives Contents monotonically toward zero with no equilibrium point.

Physically, this is wrong. As the gas cell deflates and the internal pressure drops toward ambient, the pressure differential at the valve should diminish. More importantly, once internal pressure equals ambient, air would begin flowing in through the valve to fill the void — the valve should be bi-directional.

I prototyped a fix: making CellHeight scale with the ratio of current volume to maximum volume (CellHeight *= GasVolume / MaxVolume), so DeltaPressure naturally diminishes to zero as the cell deflates. This addresses the existing FIXME on line 319 (FixMe: CellHeight should depend on current volume).

However, this breaks the weather balloon burst mechanism. The aircraft XML uses valve_coefficient = 100000 to simulate catastrophic burst via instantaneous gas release through the valve. With volume-scaled CellHeight, the venting slows to a crawl before complete emptying (as CellHeight → 0, DeltaPressure → 0), leaving a tiny residual that triggers a table lookup assertion in the aero model.

For the FGBallonet case, the issue is simpler: DeltaPressure = Pressure - AirPressure, and equilibrium is when Pressure = AirPressure. Since IdealPressure = Contents * R * T / MaxVolume, the equilibrium contents is:

Contents_eq = AirPressure * MaxVolume / (R * Temperature)

Clamping to max(Contents_eq, ...) instead of max(1.0, ...) would be physically correct.

What would you suggest as the right scope for this PR? I see two possible paths:

1. Minimal: Keep the current early-return approach (which handles the burst case correctly) and address the valve physics in a separate PR.

2. Combined: Fix the valve physics (dynamic CellHeight + ballonet equilibrium clamping) and add a separate burst/destruction mechanism that bypasses the valve entirely.

Happy to go either direction.

[DheerendraTomar]

@bcoconni @AirshipSim Would you please give me a direction so that I can proceed with it?

[bcoconni]

@DheerendraTomar, sorry for my late response. I have missed your previous message, my apologies.

Thanks for having investigated the topic in such depth. Your approach looks very interesting and is more in line with what I was expecting to address the topic properly. I will need a bit more time to digest your proposal (sorry not much free time at the moment). I should be able to give you some feedback later this week, no later than the coming week-end.

In the meantime, I would like @andgi to review your suggestion as well since he is the original author of the code.

[bcoconni]

Thanks for the proposal @DheerendraTomar.

. With volume-scaled CellHeight, the venting slows to a crawl before complete emptying (as CellHeight → 0, DeltaPressure → 0), leaving a tiny residual.

Yes, that's really the point.

What would you suggest as the right scope for this PR? I see two possible paths:

1. Minimal: Keep the current early-return approach (which handles the burst case correctly) and address the valve physics in a separate PR.
2. Combined: Fix the valve physics (dynamic CellHeight + ballonet equilibrium clamping) and add a separate burst/destruction mechanism that bypasses the valve entirely.

I would vote for option 2. Regarding the burst event, I would force the content to a small value 0.001 rather than opening the valve:

diff --git a/aircraft/weather-balloon/weather-balloon.xml b/aircraft/weather-balloon/weather-balloon.xml
index 3f0557254..9e4790351 100644
--- a/aircraft/weather-balloon/weather-balloon.xml
+++ b/aircraft/weather-balloon/weather-balloon.xml
@@ -190,7 +190,6 @@
       buoyant_forces/gas-cell/volume-ft3 ge buoyant_forces/gas-cell/max_volume-ft3
       buoyant_forces/gas-cell/burst ge 0.5
     </test>
-    <output>buoyant_forces/gas-cell/valve_open</output>
    </switch>
   </channel>
 
diff --git a/scripts/weather-balloon.xml b/scripts/weather-balloon.xml
index 295f94bab..9da4ebbdc 100644
--- a/scripts/weather-balloon.xml
+++ b/scripts/weather-balloon.xml
@@ -31,6 +31,7 @@
         <property>velocities/vt-fps</property>
         <property>metrics/radius-ft</property>
       </notify>
+      <set name="buoyant_forces/gas-cell/contents-mol" value="0.001"/>
       <condition> buoyant_forces/gas-cell/burst ge 0.5 </condition>
     </event>

With these values, the weather balloon hits the ground at t=6752.38s instead of t=6752.8s

[DheerendraTomar]

@bcoconni @andgi I think you would need to approve for the pipeline to run.

[DheerendraTomar]

@andgi can you please trigger the pipeline?

[seanmcleod70]

@DheerendraTomar I've approved the CI workflows.

[codecov]

Codecov Report

❌ Patch coverage is 0% with 6 lines in your changes missing coverage. Please review.
✅ Project coverage is 25.10%. Comparing base (0d44e12) to head (7f2d75a).

Files with missing lines	Patch %	Lines
src/models/FGGasCell.cpp	0.00%	6 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1394   +/-   ##
=======================================
  Coverage   25.10%   25.10%           
=======================================
  Files         171      171           
  Lines       18843    18843           
=======================================
  Hits         4731     4731           
  Misses      14112    14112           

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[bcoconni]

Thanks for the update and sorry of having taken so much time to review it.

You will see I have a number of comments:

• The code in the if (Contents <= 0.0) sections is redundant with the existing computation path which leads to the same values than the ones you are forcing. If this code is for optimization purpose then please remove it because:
  1. Do not optimize for edge cases. Empty gas cells seldomly occur and this PR is meant to make JSBSim behave well including for corner cases.
  2. FGGasCell is not the place where JSBSim is spending much CPU cycles anyway
  3. It makes the code harder to read.
• I have some doubts about using AirPressure to compute ContentsEquilibrium which is supposed to be a max value (and not a min value) so it would cap the content to a value that is too low when Pressure > AirPressure.

[DheerendraTomar]

Thanks for the detailed review, @bcoconni — it moved the PR to a better place. I've reworked it and the diff is now down to 3 files, +12/−10.

1. The redundant if (Contents <= 0.0) blocks — removed (both)

You were right, they were redundant. I removed both the gas-cell and the ballonet block. I confirmed — by reading the code and by running TestLighterThanAir under SIGFPE — that an empty cell flows through the normal path correctly: the temperature update is already guarded by if (Contents > 0), Pressure collapses to AirPressure, Volume to BallonetsVolume, and mass = 0 zeroes the inertia and moments — exactly the values the special case was forcing. No early return needed.

2. ContentsEquilibrium — Pressure vs AirPressure

Looking into your question, using Pressure turns out to be self-referential: in that branch Pressure = max(IdealPressure, ParentPressure), and whenever IdealPressure ≥ ParentPressure we have Pressure = Contents·R·T/MaxVolume, so Pressure·MaxVolume/(R·T) reduces back to Contents itself — i.e. the clamp would forbid any venting at all. So neither
Pressure nor AirPressure gives a clean equilibrium there. Rather than patch the value, I reverted the ballonet valving entirely to master — the questioned ContentsEquilibrium clamp is gone and FGBallonet::Calculate is now byte-identical
to master. Reasoning below.

3. Why I left the ballonet's max(1.0, …) "hack" in place

I tried the clean approach first — max(0.0, …), letting the ballonet deflate fully. While stress-testing it I found that driving a ballonet to exactly Contents == 0 produces a NaN that propagates into an aerodynamic table lookup and trips assert(Factor >= 0.0 && Factor <= 1.0) in FGTable::GetValue (FGTable.cpp:513). Notably:

• Tiny-but-nonzero contents (down to 1e-12) are fine — only exactly 0.0 fails.
• master does not hit this, because max(1.0, …) floors valved content at 1 mol so a ballonet never reaches 0 through the valve. My max(0.0, …) made 0 reachable and newly exposed the assertion — a regression I'd be introducing.

So that old FIXME: Too small values of Contents sometimes leads to NaN floor is actually guarding a real downstream fragility in the airship aero path. Letting a ballonet deflate fully (which I agree is the physically correct behaviour) needs that NaN fixed at its source first. Since this PR targets the weather-balloon burst (#654, which has no ballonets)
and the airship valving is really @andgi's domain, I've kept the ballonet untouched here. I'll open a separate issue for the empty-ballonet path with the repro.

The gas-cell valve fix (the core of the PR)

• CellHeight is scaled by GasVolume / MaxVolume, so as the cell deflates the buoyancy head → 0 and the valve self-arrests at equilibrium (Pressure → AirPressure) at a healthy, positive content — no spontaneous vacuum, addressing the 2nd-law issue you raised.
• DeltaPressure is clamped to ≥ 0 for the same reason.
• With that self-arresting valve, the max(1E-8, …) floor is no longer needed and becomes max(0.0, …).

Burst event uses 0.0, not 0.001

I initially used your 0.001, but it fails TestLighterThanAir: that test forces the burst at IC on the ground, and a tiny residual content sends the temperature update (… /(Cv·Contents·R)) unstable → negative volume → the same FGTable assertion. Setting it to exactly 0.0 takes the guarded empty-cell path and is stable. End-to-end, the balloon now
bursts at t≈4278 s and lands at t≈6751 s, matching your earlier numbers.

[bcoconni]

Looks good to me.

Note that I took the liberty to remove the line #include <algorithm> in FGGasCell.cpp since it is not needed.

[bcoconni]

The PR is now merged.

Thanks for contribution @DheerendraTomar and thanks for your patience since I took quite some time to review your proposals.