From 7e6bef4a905cbd82eac91affece6c9781d9abb8e Mon Sep 17 00:00:00 2001 From: SkyratBot <59378654+SkyratBot@users.noreply.github.com> Date: Fri, 6 May 2022 13:49:16 +0200 Subject: [PATCH] [NO GBP] Loosens gas pressure transfer calculation limit a bit. (#66553) (#13363) I initially kept it an exclusive check so funny things like division by inf doesnt get in, but I implemented more inf checks on the quadratic proc and the likes so it really shouldnt be necessary now. This should hopefully reduce cases when the target pressure hits the upper limit because the target gasmix has a miniscule amount of moles, having the result be rejected by the code, and triggering a runtime. And some minor improvements to the code so it doesnt iterate through the gaslist multiple times. Nicer code Co-authored-by: vincentiusvin <54709710+vincentiusvin@users.noreply.github.com> --- code/__DEFINES/atmospherics/atmos_core.dm | 4 +- code/__DEFINES/maths.dm | 2 +- .../atmospherics/gasmixtures/gas_mixture.dm | 40 +++++++----- code/modules/unit_tests/gas_transfer.dm | 64 ++++++++++++------- 4 files changed, 70 insertions(+), 40 deletions(-) diff --git a/code/__DEFINES/atmospherics/atmos_core.dm b/code/__DEFINES/atmospherics/atmos_core.dm index 8105920ecd7..307167cc56b 100644 --- a/code/__DEFINES/atmospherics/atmos_core.dm +++ b/code/__DEFINES/atmospherics/atmos_core.dm @@ -179,4 +179,6 @@ #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 +#define ATMOS_PRESSURE_APPROXIMATION_ITERATIONS 20 +/// We deal with big numbers and a lot of math, things are bound to get imprecise. Take this traveller. +#define ATMOS_PRESSURE_ERROR_TOLERANCE 0.01 diff --git a/code/__DEFINES/maths.dm b/code/__DEFINES/maths.dm index 6a2ae4fa531..91831182b99 100644 --- a/code/__DEFINES/maths.dm +++ b/code/__DEFINES/maths.dm @@ -1,5 +1,5 @@ // Remove these once we have Byond implementation. -#define ISNAN(a) (!(a==a)) +#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 diff --git a/code/modules/atmospherics/gasmixtures/gas_mixture.dm b/code/modules/atmospherics/gasmixtures/gas_mixture.dm index c7d96d7c721..7cf75c13296 100644 --- a/code/modules/atmospherics/gasmixtures/gas_mixture.dm +++ b/code/modules/atmospherics/gasmixtures/gas_mixture.dm @@ -578,15 +578,20 @@ GLOBAL_LIST_INIT(gaslist_cache, init_gaslist_cache()) * - 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)) + // So we dont need to iterate the gaslist multiple times. + var/our_moles = total_moles() + var/output_moles = output_air.total_moles() + var/output_pressure = output_air.return_pressure() + + if(our_moles <= 0 || temperature <= 0) return FALSE var/pressure_delta = 0 - if((output_air.temperature <= 0) || (output_air.total_moles() <= 0)) + if(output_air.temperature <= 0 || output_moles <= 0) ignore_temperature = TRUE pressure_delta = target_pressure else - pressure_delta = target_pressure - output_air.return_pressure() + pressure_delta = target_pressure - output_pressure if(pressure_delta < 0.01 || gas_pressure_minimum_transfer(output_air) > target_pressure) return FALSE @@ -596,12 +601,17 @@ GLOBAL_LIST_INIT(gaslist_cache, init_gaslist_cache()) // 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) + /// The PV/R part in the equation we will use later. Counted early because pv/(r*t) might not be equal to pv/r/t, messing our lower and upper limit. + var/pvr = pv / R_IDEAL_GAS_EQUATION // 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. + + /// This is the true lower limit, but numbers still can get lower than this due to floats. + var/lower_limit = max((pvr / max(temperature, output_air.temperature)) - output_moles, 0) + var/upper_limit = (pvr / min(temperature, output_air.temperature)) - output_moles // In theory this should never go below zero, the pressure_delta check above should account for this. + + lower_limit = max(lower_limit - ATMOS_PRESSURE_ERROR_TOLERANCE, 0) + upper_limit += ATMOS_PRESSURE_ERROR_TOLERANCE /* * We have PV=nRT as a nice formula, we can rearrange it into nT = PV/R @@ -627,17 +637,14 @@ GLOBAL_LIST_INIT(gaslist_cache, init_gaslist_cache()) // Our thermal energy and moles var/w2 = thermal_energy() - var/n2 = total_moles() + var/n2 = our_moles var/c2 = heat_capacity() // Target thermal energy and moles var/w1 = output_air.thermal_energy() - var/n1 = output_air.total_moles() + var/n1 = output_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 @@ -659,8 +666,8 @@ GLOBAL_LIST_INIT(gaslist_cache, init_gaslist_cache()) /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 + solution = max(SolveQuadratic(a, b, c)) + if(solution > lower_limit && solution < upper_limit) //SolveQuadratic can return empty lists 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 @@ -670,8 +677,9 @@ GLOBAL_LIST_INIT(gaslist_cache, init_gaslist_cache()) /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 + // We start at the extrema of the equation, added by a number. + // This way we will hopefully always converge on the positive root, while starting at a reasonable number. + solution = (-b / (2 * a)) + 200 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) diff --git a/code/modules/unit_tests/gas_transfer.dm b/code/modules/unit_tests/gas_transfer.dm index 32f973d7244..d403eafeb65 100644 --- a/code/modules/unit_tests/gas_transfer.dm +++ b/code/modules/unit_tests/gas_transfer.dm @@ -2,31 +2,51 @@ /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) + for (var/hot_test in list(1e4, 1e6, 1e8, 1e10, 1e12)) + nob_to_trit(hot_test, hot_test, 50, T20C, max(2500, hot_test/100)) + for (var/cold_test in list(1, 1e-2, MOLAR_ACCURACY)) + nob_to_trit(5000, T20C, cold_test, cold_test) + nob_to_trit(5000, T20C, 100, T20C, 1) - first_mix.volume = 200 - second_mix.volume = 200 +/** + * Proc to transfer x moles of x temp nob to x moles of x temp trit. + * + * Arguments: + * * nob_moles: Moles for the nob (origin) + * * nob_temp: Temp for the nob (origin) + * * trit_moles: Moles for the trit (target) + * * nob_temp: Temp for the nob (target) + * * additional_pressure: Optional proc, if unfilled transfer will be 10% of pressure. + */ +/datum/unit_test/atmospheric_gas_transfer/proc/nob_to_trit(nob_moles, nob_temp, trit_moles, trit_temp, additional_pressure) + var/datum/gas_mixture/first_mix = allocate(/datum/gas_mixture) + var/datum/gas_mixture/second_mix = allocate(/datum/gas_mixture) - 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 + first_mix.volume = 200 + second_mix.volume = 200 - 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 + ASSERT_GAS(/datum/gas/hypernoblium, first_mix) + ASSERT_GAS(/datum/gas/tritium, second_mix) + + first_mix.gases[/datum/gas/hypernoblium][MOLES] = nob_moles + first_mix.temperature = nob_temp + + second_mix.gases[/datum/gas/tritium][MOLES] = trit_moles + second_mix.temperature = trit_temp + + var/initial_pressure = second_mix.return_pressure() + // A fixed number would mean transfer is too small for high temps. So we make it scaled. + + if(isnull(additional_pressure)) + additional_pressure = first_mix.return_pressure() / 10 - /* 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 + /* 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)) + 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()].") + TEST_ASSERT(margin<=error_margin, "Failed to pump [nob_moles] moles of [nob_temp] K Nob to [trit_moles] moles of [trit_temp] K Trit, . Expected pressure = [initial_pressure+additional_pressure] +/- [error_margin]. Got [second_mix.return_pressure()].")