diff --git a/code/__DEFINES/atmospherics/atmos_core.dm b/code/__DEFINES/atmospherics/atmos_core.dm index c278b5b0fa8..ca906309a0f 100644 --- a/code/__DEFINES/atmospherics/atmos_core.dm +++ b/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 diff --git a/code/__DEFINES/maths.dm b/code/__DEFINES/maths.dm index c4cf7e10194..6a2ae4fa531 100644 --- a/code/__DEFINES/maths.dm +++ b/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 diff --git a/code/modules/atmospherics/gasmixtures/gas_mixture.dm b/code/modules/atmospherics/gasmixtures/gas_mixture.dm index c34ea2f469d..e0e547cea71 100644 --- a/code/modules/atmospherics/gasmixtures/gas_mixture.dm +++ b/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 diff --git a/code/modules/unit_tests/_unit_tests.dm b/code/modules/unit_tests/_unit_tests.dm index c84e2a4cf86..d5cc8a604bd 100644 --- a/code/modules/unit_tests/_unit_tests.dm +++ b/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" diff --git a/code/modules/unit_tests/gas_transfer.dm b/code/modules/unit_tests/gas_transfer.dm new file mode 100644 index 00000000000..32f973d7244 --- /dev/null +++ b/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()].")