From a07dd0b4caeeaed77d241a36bd85f218647a2814 Mon Sep 17 00:00:00 2001 From: Devin Kees <6895754+djkees@users.noreply.github.com> Date: Wed, 2 Sep 2026 10:39:14 -0600 Subject: [PATCH] Fix mixture entropy/gibbs_energy single-multi scaling mismatches (#162) entropy_multi had a stray /1.d3, gibbs_energy_multi's nj was never normalized by sum(weights), and both entropy overloads were missing /std_pressure in their pressure log term. Co-authored-by: Claude Sonnet 5 (cherry picked from commit 726abee8264684c82ab17bcb830c64c30c774d94) --- source/mixture.f90 | 16 +++++++++------- source/mixture_test.pf | 8 ++++---- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/source/mixture.f90 b/source/mixture.f90 index 2bd5a1c..09cab78 100644 --- a/source/mixture.f90 +++ b/source/mixture.f90 @@ -858,8 +858,8 @@ function mixture_calc_entropy_single(self, weights, temperature, pressure) resul s = s + nj*(sj - log(nj)) end if end do - s = s - n*log(pressure/n) - s = s * gas_constant !/ 1.d3 + s = s - n*log(pressure/std_pressure/n) + s = s * gas_constant end function @@ -892,11 +892,11 @@ function mixture_calc_entropy_multi(self, weights, temperatures, pressures) resu s = s + nj*sj else n = n + nj - s = s + nj*(sj - log(nj*pressures(j))) + s = s + nj*(sj - log(nj*pressures(j)/std_pressure)) end if end do s = s + n*log(n) - s = s * gas_constant / 1.d3 + s = s * gas_constant end function @@ -931,24 +931,26 @@ function mixture_calc_gibbs_energy_multi(self, weights, temperatures, pressures) ! Locals integer :: j - real(dp) :: gj, nj, n, nr, pr + real(dp) :: gj, nj, n, nr, pr, total_weight call check_array_len(size(weights), self%num_species, 'mixture_calc_gibbs_energy_multi weights') call check_array_len(size(temperatures), self%num_species, 'mixture_calc_gibbs_energy_multi temperatures') call check_array_len(size(pressures), self%num_species, 'mixture_calc_gibbs_energy_multi pressures') + total_weight = sum(weights) + ! Must pre-compute because cannot factor out of the ! sum when consitituent temperatures can vary. n = 0.0d0 do j = 1, self%num_species if (self%is_condensed(j)) cycle - nj = weights(j)/self%species(j)%molecular_weight + nj = weights(j)/self%species(j)%molecular_weight/total_weight n = n + nj end do g = 0.0d0 do j = 1, self%num_species - nj = weights(j)/self%species(j)%molecular_weight + nj = weights(j)/self%species(j)%molecular_weight/total_weight gj = self%species(j)%calc_gibbs_energy(temperatures(j)) if (self%is_condensed(j)) then g = g + nj*gj diff --git a/source/mixture_test.pf b/source/mixture_test.pf index 2489d26..2850aa0 100644 --- a/source/mixture_test.pf +++ b/source/mixture_test.pf @@ -329,16 +329,16 @@ contains t = temp(1) @assertRelativelyEqual(-1.19532599d+07, mix%calc_enthalpy(wtfracs, t), tol) @assertRelativelyEqual(-1.20517040d+07, mix%calc_energy(wtfracs, t), tol) - !@assertRelativelyEqual(-1.48300159E+07, mix%calc_gibbs_energy(wtfracs, t, p), tol) - !@assertRelativelyEqual( 9.46301312d+03, mix%calc_entropy(wtfracs, t, p), tol) + @assertRelativelyEqual(-1.48300159E+07, mix%calc_gibbs_energy(wtfracs, t, p), tol) + @assertRelativelyEqual( 9.46301312d+03, mix%calc_entropy(wtfracs, t, p), tol) ! Version 2: Allow per-species temperatures and pressures ! This enables fuels and oxidizers to have different initial states ! prior to being mixed and equilibrated. @assertRelativelyEqual(-1.19532599d+07, mix%calc_enthalpy(wtfracs, temp), tol) @assertRelativelyEqual(-1.20517040d+07, mix%calc_energy(wtfracs, temp), tol) - !@assertRelativelyEqual(-1.48300159E+07, mix%calc_gibbs_energy(wtfracs, temp, pres), tol) - !@assertRelativelyEqual( 9.46301312d+03, mix%calc_entropy(wtfracs, temp, pres), tol) + @assertRelativelyEqual(-1.48300159E+07, mix%calc_gibbs_energy(wtfracs, temp, pres), tol) + @assertRelativelyEqual( 9.46301312d+03, mix%calc_entropy(wtfracs, temp, pres), tol) end subroutine