Fix failing transport property cases (nasa#37)

Fix transport-property failures and legacy-behavior mismatches in CH4/O2 equilibrium/rocket cases, including an out-of-bounds crash path and interval-selection differences.

Author: markleader

source/equilibrium.f90

@@ -2,17 +2,20 @@ module equilibrium_test
     use funit
     use cea_equilibrium
     use cea_thermo
+    use cea_transport
     use cea_mixture
     use cea_units
     use cea_param, only: R=>gas_constant
 
     type(ThermoDB) :: all_thermo
+    type(TransportDB) :: all_transport
 
 contains
 
     @before
     subroutine setup_mixture()
         all_thermo = read_thermo('data/thermo.lib')
+        all_transport = read_transport('data/trans.lib')
     end subroutine
 
     @test
@@ -659,6 +662,36 @@ contains
 
     end subroutine
 
+    @test
+    subroutine test_equilibrium_transport_ch4_o2
+        type(Mixture) :: products
+        type(Mixture) :: reactants
+        type(EqSolver) :: solver
+        type(EqSolution) :: solution
+        character(:), allocatable :: product_names(:)
+        real(dp) :: h_reac, p_reac, weights(2)
+
+        reactants = Mixture(all_thermo, ['CH4', 'O2 '])
+        product_names = reactants%get_products(all_thermo)
+        products = Mixture(all_thermo, product_names)
+
+        solver = EqSolver(products, reactants, all_transport=all_transport)
+        solution = EqSolution(solver)
+
+        weights = reactants%weights_from_of([0.0d0, 1.0d0], [1.0d0, 0.0d0], 2.6d0)
+        h_reac = reactants%calc_enthalpy(weights, [298.15d0, 298.15d0])/R
+        p_reac = psi_to_bar(1000.0d0)
+
+        call solver%solve(solution, 'hp', h_reac, p_reac, weights)
+
+        @assertTrue(solution%converged)
+        @assertTrue(solution%viscosity > 0.0d0)
+        @assertTrue(solution%conductivity_fr > 0.0d0)
+        @assertTrue(solution%conductivity_eq > 0.0d0)
+        @assertTrue(solution%pr_fr > 0.0d0)
+        @assertTrue(solution%pr_eq > 0.0d0)
+    end subroutine
+
     ! @test
     ! subroutine test_inert_reactant
     !     type(Mixture) :: products
@@ -708,4 +741,4 @@ contains
 
     ! end subroutine
 
-end module
+end module

source/equilibrium_test.pf

@@ -3,17 +3,20 @@ module rocket_test
     use cea_rocket
     use cea_equilibrium
     use cea_thermo
+    use cea_transport
     use cea_mixture
     use cea_units
     use cea_param, only: R=>gas_constant, snl=>species_name_len
 
     type(ThermoDB) :: all_thermo
+    type(TransportDB) :: all_transport
 
 contains
 
     @before
     subroutine setup_mixture()
         all_thermo = read_thermo('data/thermo.lib')
+        all_transport = read_transport('data/trans.lib')
     end subroutine
 
     @test
@@ -448,4 +451,62 @@ contains
 
     end subroutine
 
+    @test
+    subroutine test_rocket_transport
+        type(Mixture) :: products
+        type(Mixture) :: reactants
+        type(RocketSolver) :: solver
+        type(RocketSolution) :: soln
+        character(:), allocatable :: product_names(:)
+        real(dp) :: hc, pc, weights(2)
+        real(dp), parameter :: tol = 1.0d-3
+
+        reactants = Mixture(all_thermo, ['CH4', 'O2 '])
+        product_names = reactants%get_products(all_thermo)
+        products = Mixture(all_thermo, product_names)
+
+        weights = reactants%weights_from_of([0.0d0, 1.0d0], [1.0d0, 0.0d0], 2.6d0)
+        hc = reactants%calc_enthalpy(weights, [298.15d0, 298.15d0])/R
+        pc = psi_to_bar(1000.0d0)
+
+        solver = RocketSolver(products, reactants, all_transport=all_transport)
+        soln = solver%solve(weights, pc, pi_p=[10.0d0], supar=[20.0d0], hc=hc)
+
+        ! Test convergence
+        @assertEqual(.true., soln%converged)
+
+        ! Chamber transport properties
+        @assertRelativelyEqual(3371.84d0,  soln%eq_soln(1)%T,               tol)
+        @assertRelativelyEqual(1.0391d0,   soln%eq_soln(1)%viscosity,       tol)
+        @assertRelativelyEqual(10.4826d0,  soln%eq_soln(1)%conductivity_eq, tol)
+        @assertRelativelyEqual(0.4635d0,   soln%eq_soln(1)%pr_eq,           tol)
+        @assertRelativelyEqual(4.1984d0,   soln%eq_soln(1)%conductivity_fr, tol)
+        @assertRelativelyEqual(0.6279d0,   soln%eq_soln(1)%pr_fr,           tol)
+
+        ! Throat transport properties
+        @assertRelativelyEqual(3148.65d0,  soln%eq_soln(2)%T,               tol)
+        @assertRelativelyEqual(0.98981d0,  soln%eq_soln(2)%viscosity,       tol)
+        @assertRelativelyEqual(8.7103d0,   soln%eq_soln(2)%conductivity_eq, tol)
+        @assertRelativelyEqual(0.4693d0,   soln%eq_soln(2)%pr_eq,           tol)
+        @assertRelativelyEqual(3.9367d0,   soln%eq_soln(2)%conductivity_fr, tol)
+        @assertRelativelyEqual(0.6324d0,   soln%eq_soln(2)%pr_fr,           tol)
+
+        ! Exit transport properties (pi_p=10, Ae/At=2.3019)
+        @assertRelativelyEqual(2457.73d0,  soln%eq_soln(3)%T,               tol)
+        @assertRelativelyEqual(0.82957d0,  soln%eq_soln(3)%viscosity,       tol)
+        @assertRelativelyEqual(4.2652d0,   soln%eq_soln(3)%conductivity_eq, tol)
+        @assertRelativelyEqual(0.5445d0,   soln%eq_soln(3)%pr_eq,           tol)
+        @assertRelativelyEqual(3.1312d0,   soln%eq_soln(3)%conductivity_fr, tol)
+        @assertRelativelyEqual(0.6423d0,   soln%eq_soln(3)%pr_fr,           tol)
+
+        ! Exit transport properties (supar=20, Ae/At=20)
+        @assertRelativelyEqual(1442.09d0,  soln%eq_soln(4)%T,               tol)
+        @assertRelativelyEqual(0.56764d0,  soln%eq_soln(4)%viscosity,       tol)
+        @assertRelativelyEqual(2.0853d0,   soln%eq_soln(4)%conductivity_eq, tol)
+        @assertRelativelyEqual(0.6351d0,   soln%eq_soln(4)%pr_eq,           tol)
+        @assertRelativelyEqual(1.9610d0,   soln%eq_soln(4)%conductivity_fr, tol)
+        @assertRelativelyEqual(0.6315d0,   soln%eq_soln(4)%pr_fr,           tol)
+
+    end subroutine
+
 end module

source/rocket_test.pf

@@ -15,13 +15,16 @@ module cea_transport
         ! Required items
         character(16) :: name
             !! Species name
-        integer :: num_intervals
-            !! Number of temperature intervals
+        integer :: num_eta_intervals
+            !! Number of temperature intervals for viscosity
+        integer :: num_lambda_intervals
+            !! Number of temperature intervals for thermal conductivity
 
         ! Coefficients
-        ! NOTE: This assumes temperature intervals are the same for  eta and lambda, which should be true
-        real(dp), allocatable :: T_fit(:, :)
-            !! Temperature intervals
+        real(dp), allocatable :: eta_T_fit(:, :)
+            !! Temperature intervals for viscosity fit
+        real(dp), allocatable :: lambda_T_fit(:, :)
+            !! Temperature intervals for thermal conductivity fit
         type(TransportFit), allocatable :: eta_coef(:)
             !! Viscosity curve fit coefficients
         type(TransportFit), allocatable :: lambda_coef(:)
@@ -83,10 +86,15 @@ elemental function st_calc_eta(self, T) result(eta)
         real(dp) :: eta
         integer :: i, idx
 
+        if (self%num_eta_intervals <= 0) then
+            eta = 0.0d0
+            return
+        end if
+
         ! Select temperature range
         idx = 1
-        do i = 1,self%num_intervals
-            if (T > self%T_fit(i, 1)) then
+        do i = 2, self%num_eta_intervals
+            if (self%eta_T_fit(i, 2) > 0.0d0 .and. T > self%eta_T_fit(i-1, 2)) then
                 idx = i
             end if
         end do
@@ -101,10 +109,15 @@ elemental function st_calc_lambda(self, T) result(lambda)
         real(dp) :: lambda
         integer :: i, idx
 
+        if (self%num_lambda_intervals <= 0) then
+            lambda = 0.0d0
+            return
+        end if
+
         ! Select temperature range
         idx = 1
-        do i = 1,self%num_intervals
-            if (T > self%T_fit(i, 1)) then
+        do i = 2, self%num_lambda_intervals
+            if (self%lambda_T_fit(i, 2) > 0.0d0 .and. T > self%lambda_T_fit(i-1, 2)) then
                 idx = i
             end if
         end do
@@ -119,10 +132,15 @@ elemental function bt_calc_eta(self, T) result(eta)
         real(dp) :: eta
         integer :: i, idx
 
+        if (self%num_intervals <= 0) then
+            eta = 0.0d0
+            return
+        end if
+
         ! Select temperature range
         idx = 1
-        do i = 1,self%num_intervals
-            if (T > self%T_fit(i, 1)) then
+        do i = 2, self%num_intervals
+            if (self%T_fit(i, 2) > 0.0d0 .and. T > self%T_fit(i-1, 2)) then
                 idx = i
             end if
         end do
@@ -145,7 +163,8 @@ function read_transport(filename) result(db)
         integer :: i, j               ! Index
         character(16) :: species(2)   ! Species names
         real(dp) :: tr_data(6, 3, 2)  ! Transport data
-        integer :: n_int              ! Temporary variable for number of intervals
+        integer :: n_int_eta          ! Number of viscosity intervals
+        integer :: n_int_lambda       ! Number of conductivity intervals
 
 
         call log_info('Reading transport database file: '//filename)
@@ -172,39 +191,39 @@ function read_transport(filename) result(db)
                 db%pure_species(db%num_pure) = species(1)
 
                 ! Oversize the curvefit arrays
-                allocate(db%pure_transport(db%num_pure)%T_fit(3, 2), &
+                allocate(db%pure_transport(db%num_pure)%eta_T_fit(3, 2), &
+                         db%pure_transport(db%num_pure)%lambda_T_fit(3, 2), &
                          db%pure_transport(db%num_pure)%eta_coef(3), &
                          db%pure_transport(db%num_pure)%lambda_coef(3))
 
                 db%pure_transport(db%num_pure)%name = species(1)
 
                 ! Get the number of valid intervals for these coefficients
-                n_int = 3
-                do j = 1,3
-                    if (abs(tr_data(1, j, 1)) < 1.0d-6) then
-                        n_int = j-1
-                        exit
-                    end if
-                end do
-                db%pure_transport(db%num_pure)%num_intervals = n_int
+                n_int_eta = count_transport_intervals(tr_data(:,:,1))
+                n_int_lambda = count_transport_intervals(tr_data(:,:,2))
+                db%pure_transport(db%num_pure)%num_eta_intervals = n_int_eta
+                db%pure_transport(db%num_pure)%num_lambda_intervals = n_int_lambda
 
                 ! Set the coefficients
-                do j = 1, n_int
-                    db%pure_transport(db%num_pure)%T_fit(j, :) = tr_data(:2, j, 1)
-
+                do j = 1, n_int_eta
+                    db%pure_transport(db%num_pure)%eta_T_fit(j, :) = tr_data(:2, j, 1)
                     db%pure_transport(db%num_pure)%eta_coef(j)%A    = tr_data(3, j, 1)
                     db%pure_transport(db%num_pure)%eta_coef(j)%B    = tr_data(4, j, 1)
                     db%pure_transport(db%num_pure)%eta_coef(j)%C    = tr_data(5, j, 1)
                     db%pure_transport(db%num_pure)%eta_coef(j)%D    = tr_data(6, j, 1)
+                end do
+                do j = 1, n_int_lambda
+                    db%pure_transport(db%num_pure)%lambda_T_fit(j, :) = tr_data(:2, j, 2)
 
                     db%pure_transport(db%num_pure)%lambda_coef(j)%A = tr_data(3, j, 2)
                     db%pure_transport(db%num_pure)%lambda_coef(j)%B = tr_data(4, j, 2)
                     db%pure_transport(db%num_pure)%lambda_coef(j)%C = tr_data(5, j, 2)
                     db%pure_transport(db%num_pure)%lambda_coef(j)%D = tr_data(6, j, 2)
                 end do
-                db%pure_transport(db%num_pure)%T_fit = db%pure_transport(db%num_pure)%T_fit(:n_int, :)
-                db%pure_transport(db%num_pure)%eta_coef = db%pure_transport(db%num_pure)%eta_coef(:n_int)
-                db%pure_transport(db%num_pure)%lambda_coef = db%pure_transport(db%num_pure)%lambda_coef(:n_int)
+                db%pure_transport(db%num_pure)%eta_T_fit = db%pure_transport(db%num_pure)%eta_T_fit(:n_int_eta, :)
+                db%pure_transport(db%num_pure)%lambda_T_fit = db%pure_transport(db%num_pure)%lambda_T_fit(:n_int_lambda, :)
+                db%pure_transport(db%num_pure)%eta_coef = db%pure_transport(db%num_pure)%eta_coef(:n_int_eta)
+                db%pure_transport(db%num_pure)%lambda_coef = db%pure_transport(db%num_pure)%lambda_coef(:n_int_lambda)
 
             else
                 ! Binary species
@@ -218,26 +237,20 @@ function read_transport(filename) result(db)
                 db%binary_transport(db%num_binary)%name = species
 
                 ! Get the number of valid intervals for these coefficients
-                n_int = 3
-                do j = 1,3
-                    if (abs(tr_data(1, j, 1)) < 1.0d-6) then
-                        n_int = j-1
-                        exit
-                    end if
-                end do
-                db%binary_transport(db%num_binary)%num_intervals = n_int
+                n_int_eta = count_transport_intervals(tr_data(:,:,1))
+                db%binary_transport(db%num_binary)%num_intervals = n_int_eta
 
                 ! Set the coefficients
-                do j = 1, n_int
+                do j = 1, n_int_eta
                     db%binary_transport(db%num_binary)%T_fit(j, :) = tr_data(:2, j, 1)
 
                     db%binary_transport(db%num_binary)%eta_coef(j)%A = tr_data(3, j, 1)
                     db%binary_transport(db%num_binary)%eta_coef(j)%B = tr_data(4, j, 1)
                     db%binary_transport(db%num_binary)%eta_coef(j)%C = tr_data(5, j, 1)
                     db%binary_transport(db%num_binary)%eta_coef(j)%D = tr_data(6, j, 1)
                 end do
-                db%binary_transport(db%num_binary)%T_fit = db%binary_transport(db%num_binary)%T_fit(:n_int, :)
-                db%binary_transport(db%num_binary)%eta_coef = db%binary_transport(db%num_binary)%eta_coef(:n_int)
+                db%binary_transport(db%num_binary)%T_fit = db%binary_transport(db%num_binary)%T_fit(:n_int_eta, :)
+                db%binary_transport(db%num_binary)%eta_coef = db%binary_transport(db%num_binary)%eta_coef(:n_int_eta)
 
             end if
 
@@ -252,6 +265,27 @@ function read_transport(filename) result(db)
 
     end function
 
+    integer function count_transport_intervals(visc_tfit, cond_tfit) result(n_int)
+        real(dp), intent(in) :: visc_tfit(6,3)
+        real(dp), intent(in), optional :: cond_tfit(6,3)
+        integer :: j
+
+        n_int = 0
+        do j = 1, 3
+            if (visc_tfit(2, j) > 0.0d0) then
+                n_int = j
+                cycle
+            end if
+            if (present(cond_tfit)) then
+                if (cond_tfit(2, j) > 0.0d0) then
+                    n_int = j
+                    cycle
+                end if
+            end if
+            exit
+        end do
+    end function
+
     function get_mixture_transport(all_transport, products, ions) result(transport_db)
         ! Get the transport database for the subset of species in the mixture
         type(TransportDB), intent(in) :: all_transport