fixed liquid density values (#102)

Author: markleader

source/bind/c/bindc.F90

@@ -1151,9 +1151,11 @@ function cea_mixture_calc_property_tp(mptr, prop_type, len_weights, weights, &
         end if
         select case(prop_type)
             case (CEA_VOLUME)
-                prop_value = 1.d-2*mix%calc_pressure(weights(:ns)/sum(weights(:ns)), temperature)/pressure
+                prop_value = 1.d-2*mix%calc_pressure( &
+                    weights(:ns)/sum(weights(:ns)), temperature, include_condensed=.true.)/pressure
             case (CEA_DENSITY)
-                prop_value = 1.d2*pressure/mix%calc_pressure(weights(:ns)/sum(weights(:ns)), temperature)
+                prop_value = 1.d2*pressure/mix%calc_pressure( &
+                    weights(:ns)/sum(weights(:ns)), temperature, include_condensed=.true.)
             case (CEA_ENTROPY)
                 prop_value = mix%calc_entropy(weights(:ns), temperature, pressure)
             case (CEA_GIBBS_ENERGY)
@@ -1186,9 +1188,11 @@ function cea_mixture_calc_property_tp_multitemp(mptr, prop_type, len_weights, we
         end if
         select case(prop_type)
             case (CEA_VOLUME)
-                prop_value = mix%calc_pressure(weights(:ns)/sum(weights(:ns)), temperatures(:ns))/pressure
+                prop_value = 1.d-2*mix%calc_pressure( &
+                    weights(:ns)/sum(weights(:ns)), temperatures(:ns), include_condensed=.true.)/pressure
             case (CEA_DENSITY)
-                prop_value = pressure/mix%calc_pressure(weights(:ns)/sum(weights(:ns)), temperatures(:ns))
+                prop_value = 1.d2*pressure/mix%calc_pressure( &
+                    weights(:ns)/sum(weights(:ns)), temperatures(:ns), include_condensed=.true.)
             case (CEA_ENTROPY)
                 block
                     real(wp) :: pressures(ns)

source/bind/python/tests/test_mixture.py

@@ -32,3 +32,47 @@ def test_calc_property_rejects_unknown_type(cea_module):
     weights = np.array([0.4, 0.6], dtype=np.float64)
     with pytest.raises(ValueError):
         mix.calc_property(999999, weights, temperature=3000.0)
+
+
+@pytest.mark.parametrize(
+    ("species", "weights", "temperatures", "expected_density"),
+    [
+        (["O2(L)"], np.array([1.0]), np.array([90.0]), 0.2907792655383313),
+        (
+            ["CH4(L)", "O2(L)"],
+            np.array([0.2, 0.8]),
+            np.array([111.0, 90.0]),
+            0.22505975471169806,
+        ),
+    ],
+)
+def test_calc_liquid_density_is_finite(
+    cea_module, species, weights, temperatures, expected_density
+):
+    mix = cea_module.Mixture(species)
+
+    density = mix.calc_property(
+        cea_module.DENSITY, weights, temperature=temperatures, pressure=68.0
+    )
+    volume = mix.calc_property(
+        cea_module.VOLUME, weights, temperature=temperatures, pressure=68.0
+    )
+
+    assert np.isfinite(density)
+    assert density > 0.0
+    assert density == pytest.approx(expected_density)
+    assert density * volume == pytest.approx(1.0)
+
+
+def test_multitemperature_density_matches_single_temperature(cea_module):
+    mix = cea_module.Mixture(["O2(L)"])
+    weights = np.array([1.0])
+
+    density_single = mix.calc_property(
+        cea_module.DENSITY, weights, temperature=90.0, pressure=68.0
+    )
+    density_multi = mix.calc_property(
+        cea_module.DENSITY, weights, temperature=np.array([90.0]), pressure=68.0
+    )
+
+    assert density_multi == pytest.approx(density_single)

source/mixture.f90

@@ -1068,52 +1068,66 @@ function mixture_calc_frozen_cv_multi(self, weights, temperatures) result(cv)
 
     end function
 
-    function mixture_calc_pressure_single(self, densities, temperature) result(p)
+    function mixture_calc_pressure_single(self, densities, temperature, include_condensed) result(p)
+        ! Condensed species normally have no ideal-gas pressure contribution. The optional
+        ! override supports fixed-mixture volume and density queries for reactants.
 
         ! Arguments
         class(Mixture), intent(in) :: self
         real(dp), intent(in) :: densities(:)
         real(dp), intent(in) :: temperature
+        logical, intent(in), optional :: include_condensed
 
         ! Return
         real(dp) :: p
 
         ! Locals
         integer :: j
         real(dp) :: n, nj
+        logical :: include_condensed_
 
         call check_array_len(size(densities), self%num_species, 'mixture_calc_pressure_single densities')
 
+        include_condensed_ = .false.
+        if (present(include_condensed)) include_condensed_ = include_condensed
+
         n = 0.0d0
         do j = 1, self%num_species
-            if (self%is_condensed(j)) cycle
+            if (self%is_condensed(j) .and. .not. include_condensed_) cycle
             nj = densities(j)/self%species(j)%molecular_weight
             n = n + nj
         end do
         p = n * gas_constant * temperature
 
     end function
 
-    function mixture_calc_pressure_multi(self, densities, temperatures) result(p)
+    function mixture_calc_pressure_multi(self, densities, temperatures, include_condensed) result(p)
+        ! Condensed species normally have no ideal-gas pressure contribution. The optional
+        ! override supports fixed-mixture volume and density queries for reactants.
 
         ! Arguments
         class(Mixture), intent(in) :: self
         real(dp), intent(in) :: densities(:)
         real(dp), intent(in) :: temperatures(:)
+        logical, intent(in), optional :: include_condensed
 
         ! Return
         real(dp) :: p
 
         ! Locals
         integer :: j
         real(dp) :: nT, nj
+        logical :: include_condensed_
 
         call check_array_len(size(densities), self%num_species, 'mixture_calc_pressure_multi densities')
         call check_array_len(size(temperatures), self%num_species, 'mixture_calc_pressure_multi temperatures')
 
+        include_condensed_ = .false.
+        if (present(include_condensed)) include_condensed_ = include_condensed
+
         nT = 0.0d0
         do j = 1, self%num_species
-            if (self%is_condensed(j)) cycle
+            if (self%is_condensed(j) .and. .not. include_condensed_) cycle
             nj = densities(j)/self%species(j)%molecular_weight
             nT = nT + nj * temperatures(j)
         end do

source/mixture_test.pf

@@ -292,6 +292,21 @@ contains
         @assertRelativelyEqual(0.39655431d3, h0/gas_constant, tol)
     end subroutine
 
+    @test
+    subroutine test_calc_pressure_including_condensed
+        type(Mixture) :: mix
+        real(dp), parameter :: tol = 1.d-12
+        real(dp) :: actual, expected
+
+        mix = Mixture(all_thermo, ['CH4(L)', 'O2(L) '])
+
+        expected = gas_constant*(0.2d0*111.0d0/mix%species(1)%molecular_weight + &
+                                 0.8d0*90.0d0/mix%species(2)%molecular_weight)
+        @assertEqual(0.0d0, mix%calc_pressure([0.2d0, 0.8d0], [111.0d0, 90.0d0]))
+        actual = mix%calc_pressure([0.2d0, 0.8d0], [111.0d0, 90.0d0], include_condensed=.true.)
+        @assertRelativelyEqual(expected, actual, tol)
+    end subroutine
+
     @test
     subroutine test_calc_properties
         ! Taken from output of RP-1311 Example Case 14