[MIRROR] Pump pressure calculation now respects temperature. [MDB IGNORE] (#10812)

* Pump pressure calculation now respects temperature. (#62320)

Ever pump hot nob into a tank of oxygen, see the pressure rise to 8000 kpa for a few seconds and have that deafening explosion sound play right after? This attempts to eliminate those by making the pressure transfer proc thingy in gas_mixture account for resulting temperature. The old one is preserved if the temperature difference isn't too big because it's most definitely more performant. The minimum temperature difference of 5 is completely arbitrary and I won't mind changing if asked. Math is a bit messy by necessity but I hope the comment is adequate.

While I'm at it i also changed the return value of those procs to return a removed gasmix instead of directly merging. I sort of needed this for a clean implementation of assume_air in #62284. Will need more testing though. Found a better solution.

Writes a small unit test to check for this too.

Cuts down on one unintended cause of explosions.

* Pump pressure calculation now respects temperature.

Co-authored-by: vincentiusvin <54709710+vincentiusvin@users.noreply.github.com>

code/__DEFINES/atmospherics/atmos_core.dm

@@ -193,3 +193,6 @@
 #define MAKE_ACTIVE 2
 ///Disable excited group
 #define KILL_EXCITED 3
+
+/// How many maximum iterations do we allow the Newton-Raphson approximation for gas pressure to do.
+#define ATMOS_PRESSURE_APPROXIMATION_ITERATIONS 10

code/__DEFINES/maths.dm

@@ -1,3 +1,9 @@
+// Remove these once we have Byond implementation.
+#define ISNAN(a) (!(a==a))
+#define ISINF(a) (!ISNAN(a) && ISNAN(a-a))
+#define IS_INF_OR_NAN(a) (ISNAN(a-a))
+// Aight dont remove the rest
+
 // Credits to Nickr5 for the useful procs I've taken from his library resource.
 // This file is quadruple wrapped for your pleasure
 // (
@@ -95,7 +101,7 @@
 	. = list()
 	var/d = b*b - 4 * a * c
 	var/bottom  = 2 * a
-	if(d < 0)
+	if(d < 0 || IS_INF_OR_NAN(d) || IS_INF_OR_NAN(bottom))
 		return
 	var/root = sqrt(d)
 	. += (-b + root) / bottom

code/modules/atmospherics/gasmixtures/gas_mixture.dm

@@ -534,11 +534,19 @@ GLOBAL_LIST_INIT(gaslist_cache, init_gaslist_cache())
 	if(.) //If we changed the mix to any degree
 		garbage_collect()
 
-///Takes the amount of the gas you want to PP as an argument
-///So I don't have to do some hacky switches/defines/magic strings
-///eg:
-///Plas_PP = get_partial_pressure(gas_mixture.plasma)
-///O2_PP = get_partial_pressure(gas_mixture.oxygen)
+
+/**
+ * Takes the amount of the gas you want to PP as an argument
+ * So I don't have to do some hacky switches/defines/magic strings
+ * eg:
+ * Plas_PP = get_partial_pressure(gas_mixture.plasma)
+ * O2_PP = get_partial_pressure(gas_mixture.oxygen)
+ * get_breath_partial_pressure(gas_pp) --> gas_pp/total_moles()*breath_pp = pp
+ * get_true_breath_pressure(pp) --> gas_pp = pp/breath_pp*total_moles()
+ *
+ * 10/20*5 = 2.5
+ * 10 = 2.5/5*20
+ */
 
 /datum/gas_mixture/proc/get_breath_partial_pressure(gas_pressure)
 	return (gas_pressure * R_IDEAL_GAS_EQUATION * temperature) / BREATH_VOLUME
@@ -546,58 +554,168 @@ GLOBAL_LIST_INIT(gaslist_cache, init_gaslist_cache())
 /datum/gas_mixture/proc/get_true_breath_pressure(partial_pressure)
 	return (partial_pressure * BREATH_VOLUME) / (R_IDEAL_GAS_EQUATION * temperature)
 
-///Mathematical proofs:
-/**
-get_breath_partial_pressure(gas_pp) --> gas_pp/total_moles()*breath_pp = pp
-get_true_breath_pressure(pp) --> gas_pp = pp/breath_pp*total_moles()
+/** 
+ * Counts how much pressure will there be if we impart MOLAR_ACCURACY amounts of our gas to the output gasmix. 
+ * We do all of this without actually transferring it so dont worry about it changing the gasmix.
+ * Returns: Resulting pressure (number).
+ * Args: 
+ * - output_air (gasmix).
+ */ 
+/datum/gas_mixture/proc/gas_pressure_minimum_transfer(datum/gas_mixture/output_air)
+	var/resulting_energy = output_air.thermal_energy() + (MOLAR_ACCURACY / total_moles() * thermal_energy())
+	var/resulting_capacity = output_air.heat_capacity() + (MOLAR_ACCURACY / total_moles() * heat_capacity())
+	return (output_air.total_moles() + MOLAR_ACCURACY) * R_IDEAL_GAS_EQUATION * (resulting_energy / resulting_capacity) / output_air.volume
 
-10/20*5 = 2.5
-10 = 2.5/5*20
-**/
+
+/** Returns the amount of gas to be pumped to a specific container.
+ * Args:
+ * - output_air. The gas mix we want to pump to.
+ * - target_pressure. The target pressure we want.
+ * - ignore_temperature. Returns a cheaper form of gas calculation, useful if the temperature difference between the two gasmixes is low or nonexistant.
+ */
+/datum/gas_mixture/proc/gas_pressure_calculate(datum/gas_mixture/output_air, target_pressure, ignore_temperature = FALSE)
+	if((total_moles() <= 0) || (temperature <= 0))
+		return FALSE
+
+	var/pressure_delta = 0
+	if((output_air.temperature <= 0) || (output_air.total_moles() <= 0))
+		ignore_temperature = TRUE
+		pressure_delta = target_pressure
+	else
+		pressure_delta = target_pressure - output_air.return_pressure()
+
+	if(pressure_delta < 0.01 || gas_pressure_minimum_transfer(output_air) > target_pressure)
+		return FALSE
+
+	if(ignore_temperature)
+		return (pressure_delta*output_air.volume)/(temperature * R_IDEAL_GAS_EQUATION)
+
+	// Lower and upper bound for the moles we must transfer to reach the pressure. The answer is bound to be here somewhere.
+	var/pv = target_pressure * output_air.volume
+	var/rt_low = R_IDEAL_GAS_EQUATION * max(temperature, output_air.temperature) // Low refers to the resulting mole, this number is actually higher.
+	var/rt_high = R_IDEAL_GAS_EQUATION * min(temperature, output_air.temperature)
+	// These works by assuming our gas has extremely high heat capacity
+	// and the resultant gasmix will hit either the highest or lowest temperature possible.
+	var/lower_limit = max((pv / rt_low) - output_air.total_moles(), 0)
+	var/upper_limit = (pv / rt_high) - output_air.total_moles() // In theory this should never go below zero, the pressure_delta check above should account for this.
+
+	/*
+	 * We have PV=nRT as a nice formula, we can rearrange it into nT = PV/R
+	 * But now both n and T can change, since any incoming moles also change our temperature.
+	 * So we need to unify both our n and T, somehow.
+	 * 
+	 * We can rewrite T as (our old thermal energy + incoming thermal energy) divided by (our old heat capacity + incoming heat capacity)
+	 * T = (W1 + n/N2 * W2) / (C1 + n/N2 * C2). C being heat capacity, W being work, N being total moles.
+	 * 
+	 * In total we now have our equation be: (N1 + n) * (W1 + n/N2 * W2) / (C1 + n/N2 * C2) = PV/R
+	 * Now you can rearrange this and find out that it's a quadratic equation and pretty much solvable with the formula. Will be a bit messy though.
+	 * 
+	 * W2/N2n^2 + 
+	 * (N1*W2/N2)n + W1n - ((PV/R)*C2/N2)n + 
+	 * (-(PV/R)*C1) + N1W1 = 0
+	 * 
+	 * We will represent each of these terms with A, B, and C. A for the n^2 part, B for the n^1 part, and C for the n^0 part.
+	 * We then put this into the famous (-b +/- sqrt(b^2-4ac)) / 2a formula.
+	 * 
+	 * Oh, and one more thing. By "our" we mean the gasmix in the argument. We are the incoming one here. We are number 2, target is number 1.
+	 * If all this counting fucks up, we revert first to Newton's approximation, then the old simple formula.
+	 */
+
+	// Our thermal energy and moles
+	var/w2 = thermal_energy()
+	var/n2 = total_moles()
+	var/c2 = heat_capacity()
+	
+	// Target thermal energy and moles
+	var/w1 = output_air.thermal_energy()
+	var/n1 = output_air.total_moles()
+	var/c1 = output_air.heat_capacity()
+
+	/// The PV/R part in our equation.
+	var/pvr = pv / R_IDEAL_GAS_EQUATION
+	
+	/// x^2 in the quadratic
+	var/a_value = w2/n2
+	/// x^1 in the quadratic
+	var/b_value = ((n1*w2)/n2) + w1 - (pvr*c2/n2)
+	/// x^0 in the quadratic
+	var/c_value = (-1*pvr*c1) + n1 * w1
+	
+	. = gas_pressure_quadratic(a_value, b_value, c_value, lower_limit, upper_limit)
+	if(.)
+		return
+	. = gas_pressure_approximate(a_value, b_value, c_value, lower_limit, upper_limit)
+	if(.)
+		return
+	// Inaccurate and will probably explode but whatever.
+	return (pressure_delta*output_air.volume)/(temperature * R_IDEAL_GAS_EQUATION)
+
+/// Actually tries to solve the quadratic equation.
+/// Do mind that the numbers can get very big and might hit BYOND's single point float limit.
+/datum/gas_mixture/proc/gas_pressure_quadratic(a, b, c, lower_limit, upper_limit)
+	var/solution
+	if(!IS_INF_OR_NAN(a) && !IS_INF_OR_NAN(b) && !IS_INF_OR_NAN(c))
+		solution = max(SolveQuadratic(a, b, c)) 
+		if((solution > lower_limit) && (solution < upper_limit)) //SolveQuadratic can return nulls so be careful here
+			return solution
+	stack_trace("Failed to solve pressure quadratic equation. A: [a]. B: [b]. C:[c]. Current value = [solution]. Expected lower limit: [lower_limit]. Expected upper limit: [upper_limit].")
+	return FALSE
+
+/// Approximation of the quadratic equation using Newton-Raphson's Method.
+/// We use the slope of an approximate value to get closer to the root of a given equation.
+/datum/gas_mixture/proc/gas_pressure_approximate(a, b, c, lower_limit, upper_limit)
+	var/solution
+	if(!IS_INF_OR_NAN(a) && !IS_INF_OR_NAN(b) && !IS_INF_OR_NAN(c))
+		/// We need to start off at a reasonably good estimate. For very big numbers the amount of moles is most likely small so better start with lower_limit.
+		solution = lower_limit
+		for (var/iteration in 1 to ATMOS_PRESSURE_APPROXIMATION_ITERATIONS)
+			var/diff = (a*solution**2 + b*solution + c) / (2*a*solution + b) // f(sol) / f'(sol)
+			solution -= diff // xn+1 = xn - f(sol) / f'(sol)
+			if(abs(diff) < MOLAR_ACCURACY && (solution > lower_limit) && (solution < upper_limit))
+				return solution
+	stack_trace("Newton's Approximation for pressure failed after [ATMOS_PRESSURE_APPROXIMATION_ITERATIONS] iterations. A: [a]. B: [b]. C:[c]. Current value: [solution]. Expected lower limit: [lower_limit]. Expected upper limit: [upper_limit].")
+	return FALSE
 
 /// Pumps gas from src to output_air. Amount depends on target_pressure
 /datum/gas_mixture/proc/pump_gas_to(datum/gas_mixture/output_air, target_pressure, specific_gas = null)
-	var/output_starting_pressure = output_air.return_pressure()
+	var/temperature_delta = abs(temperature - output_air.temperature)
+	var/datum/gas_mixture/removed
+	var/transfer_moles
 
-	if((target_pressure - output_starting_pressure) < 0.01)
-		//No need to pump gas if target is already reached!
+	if(specific_gas)
+		// This is necessary because the specific heat capacity of a gas might be different from our gasmix.
+		var/datum/gas_mixture/temporary = remove_specific_ratio(specific_gas, 1)
+		transfer_moles = temporary.gas_pressure_calculate(output_air, target_pressure, temperature_delta <= 5)
+		removed = temporary.remove_specific(specific_gas, transfer_moles)
+	else
+		transfer_moles = gas_pressure_calculate(output_air, target_pressure, temperature_delta <= 5)
+		removed = remove(transfer_moles)
+	
+	if(!removed)
 		return FALSE
 
-	//Calculate necessary moles to transfer using PV=nRT
-	if((total_moles() > 0) && (temperature>0))
-		var/pressure_delta = target_pressure - output_starting_pressure
-		var/transfer_moles = (pressure_delta*output_air.volume)/(temperature * R_IDEAL_GAS_EQUATION)
-
-		//Actually transfer the gas
-		if(specific_gas)
-			var/datum/gas_mixture/removed = remove_specific(specific_gas, transfer_moles)
-			output_air.merge(removed)
-			return TRUE
-		var/datum/gas_mixture/removed = remove(transfer_moles)
-		output_air.merge(removed)
-		return TRUE
-	return FALSE
+	output_air.merge(removed)
+	return TRUE
 
 /// Releases gas from src to output air. This means that it can not transfer air to gas mixture with higher pressure.
 /datum/gas_mixture/proc/release_gas_to(datum/gas_mixture/output_air, target_pressure, rate=1)
 	var/output_starting_pressure = output_air.return_pressure()
 	var/input_starting_pressure = return_pressure()
-
+	
+	//Need at least 10 KPa difference to overcome friction in the mechanism
 	if(output_starting_pressure >= min(target_pressure,input_starting_pressure-10))
-		//No need to pump gas if target is already reached or input pressure is too low
-		//Need at least 10 KPa difference to overcome friction in the mechanism
+		return FALSE
+	//Can not have a pressure delta that would cause output_pressure > input_pressure
+	target_pressure = output_starting_pressure + min(target_pressure - output_starting_pressure, (input_starting_pressure - output_starting_pressure)/2)
+	var/temperature_delta = abs(temperature - output_air.temperature)
+	
+	var/transfer_moles = gas_pressure_calculate(output_air, target_pressure, temperature_delta <= 5)
+
+	//Actually transfer the gas
+	var/datum/gas_mixture/removed = remove(transfer_moles * rate)
+	
+	if(!removed)
 		return FALSE
 
-	//Calculate necessary moles to transfer using PV = nRT
-	if((total_moles() > 0) && (temperature>0))
-		var/pressure_delta = min(target_pressure - output_starting_pressure, (input_starting_pressure - output_starting_pressure)/2)
-		//Can not have a pressure delta that would cause output_pressure > input_pressure
-
-		var/transfer_moles = (pressure_delta*output_air.volume)/(temperature * R_IDEAL_GAS_EQUATION)
-
-		//Actually transfer the gas
-		var/datum/gas_mixture/removed = remove(transfer_moles * rate)
-		output_air.merge(removed)
-
-		return TRUE
-	return FALSE
+	output_air.merge(removed)
+	return TRUE

code/modules/unit_tests/_unit_tests.dm

@@ -67,6 +67,7 @@
 #include "egg_glands.dm"
 #include "emoting.dm"
 #include "food_edibility_check.dm"
+#include "gas_transfer.dm"
 #include "greyscale_config.dm"
 #include "heretic_knowledge.dm"
 #include "holidays.dm"

code/modules/unit_tests/gas_transfer.dm

@@ -0,0 +1,32 @@
+/// Test to make sure the pressure pumping proc used by things like portable pumps, pressure pumps, etc actually work.
+/datum/unit_test/atmospheric_gas_transfer
+
+/datum/unit_test/atmospheric_gas_transfer/Run()
+	for (var/tempNmoles in list(1e4, 1e6, 1e8, 1e10, 1e12))
+		var/datum/gas_mixture/first_mix = allocate(/datum/gas_mixture)
+		var/datum/gas_mixture/second_mix = allocate(/datum/gas_mixture)
+
+		first_mix.volume = 200
+		second_mix.volume = 200
+
+		ASSERT_GAS(/datum/gas/hypernoblium, first_mix)
+		ASSERT_GAS(/datum/gas/tritium, second_mix)
+		first_mix.gases[/datum/gas/hypernoblium][MOLES] = tempNmoles
+		second_mix.gases[/datum/gas/tritium][MOLES] = 200
+		first_mix.temperature = tempNmoles
+		second_mix.temperature = T20C
+
+		var/initial_pressure = second_mix.return_pressure()
+		// A constant value would be nicer but there will be cases when even MOLAR_ACCURACY amounts would far exceed the pressure so we need to scale it somewhat.
+		var/additional_pressure = (tempNmoles / 1000) + 500
+		
+		/* ERROR MARGIN CALCULATION
+		 * We calculate how much would the pressure change if MOLAR_ACCURACY amount of hothotgas is imparted on the cold mix.
+		 * This number gets really big for very high temperatures so it's somewhat meaningless, but our main goal is to ensure the code doesn't break.
+		 */ 
+		var/error_margin = first_mix.gas_pressure_minimum_transfer(second_mix) - initial_pressure
+		
+		first_mix.pump_gas_to(second_mix, (initial_pressure + additional_pressure))
+		var/margin = abs(second_mix.return_pressure() - (initial_pressure+additional_pressure))
+
+		TEST_ASSERT(margin<=error_margin, "Gas pressure pumping test failed for [tempNmoles]. Expected pressure = [initial_pressure+additional_pressure] +/- [error_margin]. Got [second_mix.return_pressure()].")