mixture_test.pf

module mixture_test
    use funit
    use cea_param
    use cea_thermo
    use cea_mixture
    use fb_utils

    implicit none

    type(ThermoDB) :: all_thermo

contains

    @before
    subroutine setup_thermo()
        all_thermo = read_thermo('data/thermo.lib')
    end subroutine

    @test
    subroutine test_mixture_init

        type(Mixture) :: mix
        character(snl), allocatable :: species(:)
        character(enl), allocatable :: elements(:)

        allocate(species(0), elements(0))

        species  = ['H  ', 'H2 ', 'H2O', 'O  ', 'O2 ', 'OH ']
        elements = ['H', 'O']
        mix = Mixture(all_thermo, species)

        ! Check basic meta data
        @assertEqual(6, mix%num_species)
        @assertEqual(2, mix%num_elements)

        ! Verify no re-ordering of species or elements
        @assertEqual(species(3),  mix%species_names(3))
        @assertEqual(species(6),  mix%species_names(6))
        @assertEqual(elements(1), mix%element_names(1))
        @assertEqual(elements(2), mix%element_names(2))

        ! Check correct data is loaded from the ThermoDB
        @assertEqual('H',   mix%species(1)%name)
        @assertEqual('H2O', mix%species(3)%name)
        @assertEqual('O2',  mix%species(5)%name)

        ! Check species stoichiometry
        @assertEqual([2.0d0, 0.0d0], mix%stoich_matrix(2,:))
        @assertEqual([0.0d0, 2.0d0], mix%stoich_matrix(5,:))
        @assertEqual([1.0d0, 1.0d0], mix%stoich_matrix(6,:))

    end subroutine

    @test
    subroutine test_mixture_elements_specified

        type(Mixture) :: mix1, mix2

        mix1 = Mixture(all_thermo, ['N2 ', 'N2+', 'N  ',  'N+ ', 'e- '])
        mix2 = Mixture(all_thermo, ['N2 ', 'N2+', 'N  ',  'N+ ', 'e- '], ['H ', 'N ', 'O ', 'e-'])

        @assertEqual(2, mix1%num_elements)
        @assertEqual(4, mix2%num_elements)

        @assertEqual('E', mix1%element_names(1))
        @assertEqual('N', mix1%element_names(2))
        @assertEqual('N', mix2%element_names(2))
        @assertEqual('E', mix2%element_names(4))

        @assertEqual(zeros(5),                mix2%stoich_matrix(:,1))
        @assertEqual(mix1%stoich_matrix(:,2), mix2%stoich_matrix(:,2))
        @assertEqual(zeros(5),                mix2%stoich_matrix(:,3))
        @assertEqual(mix1%stoich_matrix(:,1), mix2%stoich_matrix(:,4))

    end subroutine

    @test
    subroutine test_get_valence_ignores_electron
        type(Mixture) :: mix
        real(dp) :: vm(2), vp(2)
        real(dp), parameter :: tol = 1.d-14

        mix = Mixture(all_thermo, ['N2 ', 'N2+', 'N  ', 'N+ ', 'e- '], ions=.true.)
        call mix%get_valence(vm, vp)

        @assertEqual('N', mix%element_names(1))
        @assertEqual('E', mix%element_names(2))
        @assertTrue(abs(vm(1)) <= tol)
        @assertTrue(abs(vp(1)) <= tol)
        @assertTrue(abs(vm(2)) <= tol)
        @assertTrue(abs(vp(2)) <= tol)
    end subroutine

    @test
    subroutine test_sort_condensed

        type(Mixture) :: mix

        mix = Mixture(all_thermo, [character(7) :: 'H2', 'H2O', 'H2O(cr)', 'H2O(L)', 'O2 '])
        @assertEqual(5, mix%num_species)
        @assertEqual(3, mix%num_gas)
        @assertEqual(2, mix%num_condensed)
        @assertEqual([.false., .false., .true., .true., .false.], mix%is_condensed)

        mix = Mixture(all_thermo, [character(7) :: 'H2', 'H2O', 'H2O(cr)', 'H2O(L)', 'O2 '], sort_condensed=.true.)
        @assertEqual(5, mix%num_species)
        @assertEqual(3, mix%num_gas)
        @assertEqual(2, mix%num_condensed)
        @assertEqual([.false., .false., .false., .true., .true.], mix%is_condensed)
        @assertEqual('O2', mix%species(3)%name)
        @assertEqual('O2', mix%species_names(3))
        @assertEqual('H2O(L)', mix%species(5)%name)
        @assertEqual('H2O(L)', mix%species_names(5))

    end subroutine

    @test
    subroutine test_elements_from_species
        ! Verify elements_from_species operates correctly.

        type(Mixture) :: mix1, mix2
        real(dp) :: elements(3)

        mix1 = Mixture(all_thermo, ['H2', 'H ', 'O2', 'O '],       ['H', 'N', 'O'])
        mix2 = Mixture(all_thermo, ['N2', 'N ', 'O2', 'O ', 'NO'], ['H', 'N', 'O'])

        elements = mix1%elements_from_species([0.10d0, 0.00d0, 0.10d0, 0.00d0]) &
                 + mix2%elements_from_species([0.79d0, 0.00d0, 0.21d0, 0.00d0, 0.00d0])

        @assertEqual([0.20d0, 1.58d0, 0.62d0], elements)

    end subroutine

    @test
    subroutine test_mixture_thermo
        type(MixtureThermo) :: thermo

        thermo = MixtureThermo(12)
        @assertEqual(12, thermo%num_species)
        @assertEqual(12, size(thermo%cp))
        @assertEqual(12, size(thermo%enthalpy))
        @assertEqual(12, size(thermo%entropy))

        @assertEqual(empty_dp, thermo%cp(10))
        @assertEqual(empty_dp, thermo%enthalpy(7))
        @assertEqual(empty_dp, thermo%entropy(3))

    end subroutine

    @test
    subroutine test_calc_thermo
        type(Mixture) :: mix
        type(MixtureThermo) :: thermo
        real(dp), parameter :: T = 1200.0d0
        mix = Mixture(all_thermo, ['N2', 'O2'])
        call mix%calc_thermo(thermo, T)
        @assertEqual(2, thermo%num_species)
        @assertEqual(mix%species(1)%calc_cp(T), thermo%cp(1))
        @assertEqual(mix%species(2)%calc_enthalpy(T), T*thermo%enthalpy(2))
    end subroutine

    @test
    subroutine test_get_products
        type(Mixture) :: mix
        character(snl), allocatable :: product_names(:)
        allocate(product_names(0))
        mix = Mixture(all_thermo, ['N2', 'O2'])
        product_names = mix%get_products(all_thermo)
        @assertEqual(13, size(product_names))
        @assertEqual('N',    product_names(1))
        @assertEqual('N2',   product_names(2))
        @assertEqual('N2O',  product_names(3))
        @assertEqual('N2O3', product_names(4))
        @assertEqual('N2O4', product_names(5))
        @assertEqual('N2O5', product_names(6))
        @assertEqual('N3',   product_names(7))
        @assertEqual('NO',   product_names(8))
        @assertEqual('NO2',  product_names(9))
        @assertEqual('NO3',  product_names(10))
        @assertEqual('O',    product_names(11))
        @assertEqual('O2',   product_names(12))
        @assertEqual('O3',   product_names(13))
    end subroutine

    @test
    subroutine test_get_products_omit
        type(Mixture) :: mix
        character(snl), allocatable :: product_names(:)
        allocate(product_names(0))
        mix = Mixture(all_thermo, ['N2', 'O2'])
        product_names = mix%get_products(all_thermo, omit=['N2O3','N2O4','N2O5','N3  ','NO2 ','NO3 ','O3  '])
        @assertEqual(6, size(product_names))
        @assertEqual('N',   product_names(1))
        @assertEqual('N2',  product_names(2))
        @assertEqual('N2O', product_names(3))
        @assertEqual('NO',  product_names(4))
        @assertEqual('O',   product_names(5))
        @assertEqual('O2',  product_names(6))
    end subroutine

    @test
    subroutine test_weights_from_of
        type(Mixture) :: mix
        real(dp) :: weights(2)
        real(dp), parameter :: tol = 1.d-6

        mix = Mixture(all_thermo, ['H2', 'O2'])
        weights = mix%weights_from_of([0.0d0, 1.0d0], [1.0d0, 0.0d0], 9.0d0)
        @assertRelativelyEqual(3.40146799999999990d0, weights(1), tol)
        @assertRelativelyEqual(30.613211999999997d0, weights(2), tol)
    end subroutine

    @test
    subroutine test_chem_eq_ratio_to_of_ratio
        type(Mixture) :: mix
        real(dp) :: of_ratio
        real(dp) :: fuel_weights(2)
        real(dp) :: oxidant_weights(2)
        real(dp), parameter :: tol = 1.d-4

        mix = Mixture(all_thermo, ['H2', 'O2'])
        of_ratio = mix%of_from_equivalence([0.0d0, 1.0d0], [1.0d0, 0.0d0], 0.5d0)
        @assertRelativelyEqual(15.87337d0, of_ratio, tol)

        ! From Example 1
        mix = Mixture(all_thermo, ['H2 ', 'Air'])
        oxidant_weights = mix%weights_from_moles([0.0d0, 1.0d0])
        fuel_weights    = mix%weights_from_moles([1.0d0, 0.0d0])
        @assertRelativelyEqual(34.29623d0, mix%of_from_equivalence(oxidant_weights, fuel_weights, 1.0d0), tol)
        @assertRelativelyEqual(22.85255d0, mix%of_from_equivalence(oxidant_weights, fuel_weights, 1.5d0), tol)

    end subroutine

    @test
    subroutine test_weight_eq_ratio_to_of_ratio
        type(Mixture) :: mix
        real(dp) :: of_ratio
        real(dp), parameter :: tol = 1.d-6

        mix = Mixture(all_thermo, ['H2', 'O2'])
        of_ratio = mix%of_from_phi([0.0d0, 1.0d0], [1.0d0, 0.0d0], 0.5d0)
        @assertRelativelyEqual(15.87337d0, of_ratio, tol)
    end subroutine

    @test
    subroutine test_moles_from_weights
        type(Mixture) :: mix
        real(dp) :: moles(2)
        real(dp), parameter :: tol = 1.d-6

        mix = Mixture(all_thermo, ['H2', 'O2'])
        moles = mix%moles_from_weights([2.01588d0, 31.99880d0])
        @assertRelativelyEqual(1.0d0, moles(1), tol)
        @assertRelativelyEqual(1.0d0, moles(2), tol)
    end subroutine

    @test
    subroutine test_weights_from_moles
        type(Mixture) :: mix
        real(dp) :: weights(2)
        real(dp), parameter :: tol = 1.d-6

        mix = Mixture(all_thermo, ['H2', 'O2'])
        weights = mix%weights_from_moles([1.0d0, 1.0d0])
        @assertRelativelyEqual(2.01588d0, weights(1), tol)
        @assertRelativelyEqual(31.99880d0, weights(2), tol)
    end subroutine

    @test
    subroutine test_get_element_amounts
        type(Mixture) :: mix
        real(dp) :: b0(2), weights(2)
        real(dp), parameter :: tol = 1.d-6

        mix = Mixture(all_thermo, ['H2', 'O2'])
        weights = mix%weights_from_moles([1.0d0, 1.0d0])
        b0 = mix%element_amounts_from_weights(weights)
        @assertRelativelyEqual(0.58798142d-1, b0(1), tol)
        @assertRelativelyEqual(0.58798142d-1, b0(2), tol)
    end subroutine

    @test
    subroutine test_calc_enthalpy
        type(Mixture) :: mix
        real(dp) :: weights(2)
        real(dp) :: h0
        real(dp), parameter :: tol = 1.d-6

        mix = Mixture(all_thermo, ['H2', 'O2'])
        weights = mix%weights_from_moles([1.0d0, 1.0d0])
        h0 = mix%calc_enthalpy(weights, [2000.0d0, 2000.0d0])
        @assertRelativelyEqual(0.39655431d3, h0/gas_constant, tol)
    end subroutine

    @test
    subroutine test_calc_properties
        ! Taken from output of RP-1311 Example Case 14

        type(Mixture) :: mix
        character(*), parameter :: species(3) = [ 'H2O   ',       'O2    ',       'H2O(L)'       ]
        real(dp), parameter     :: moles(3)   = [ 6.60251584d-01, 9.09090909d-02, 2.48839325d-01 ]
        real(dp), parameter     :: pres(3)    = [ 5.06625000d+03, 5.06625000d+03, 5.06625000d+03 ]
        real(dp), parameter     :: temp(3)    = [ 304.00d0,       304.00d0,       304.00d0       ]
        real(dp), parameter     :: tol = 1d-8
        real(dp), allocatable   :: wtfracs(:)
        real(dp) :: p,t

        mix = Mixture(all_thermo, species)
        wtfracs = mix%weights_from_moles(moles)
        wtfracs = wtfracs / sum(wtfracs)

        ! Version 1: All constituents at same temp/pres
        p = pres(1)
        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)

        ! 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)

    end subroutine

end module