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()].")