diff --git a/CHANGELOG.md b/CHANGELOG.md index 8209dd8..276a8c4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,8 +6,12 @@ All notable user-visible changes to this project are documented here. ### Changed - Legacy CLI `.out` files now echo each problem’s original input block before the corresponding equilibrium, rocket, shock, or detonation output section. +- Frozen rocket failure handling was aligned with the established partial-output behavior: successful upstream stations are retained, output truncates at the last valid station, and the existing warning text is emitted after the partial report. +- Frozen rocket stop checks now use the active thermo-fit interval for condensed species plus the existing lower-temperature frozen-fit guard, and frozen retained points now populate transport properties before output. ### Fixed +- Non-converging frozen rocket cases no longer produce empty output, overrun into invalid exit columns, or crash while writing partial results. +- Difficult frozen rocket chamber solves can now retain a usable reduced-component state after repeated singular matrices, restoring partial output for reproduced failure cases. ### Added diff --git a/source/equilibrium.f90 b/source/equilibrium.f90 index a0309b8..5b496fd 100644 --- a/source/equilibrium.f90 +++ b/source/equilibrium.f90 @@ -91,6 +91,7 @@ module cea_equilibrium procedure :: test_condensed => EqSolver_test_condensed procedure :: correct_singular => EqSolver_correct_singular procedure :: update_transport_basis => EqSolver_update_transport_basis + procedure :: compute_transport_state => EqSolver_compute_transport_state procedure :: assemble_matrix => EqSolver_assemble_matrix procedure :: post_process => EqSolver_post_process procedure :: solve => EqSolver_solve @@ -739,6 +740,15 @@ subroutine EqSolver_restore_reduced_elements(self, soln, num_swaps, swap_from, s self%reduced_elements = 0 end subroutine + subroutine EqSolver_compute_transport_state(self, soln) + class(EqSolver), intent(in) :: self + type(EqSolution), intent(inout) :: soln + + if (.not. self%transport) return + call self%update_transport_basis(soln) + call compute_transport_properties(self, soln) + end subroutine + function EqSolver_compute_damped_update_factor(self, soln) result(lambda) ! Compute the damped update factor, lambda, for the Newton solver @@ -1865,7 +1875,7 @@ subroutine EqSolver_test_condensed(self, soln, iter, singular_index) end subroutine - subroutine EqSolver_correct_singular(self, soln, iter, ierr, singular_index, reduced_from, reduced_to) + subroutine EqSolver_correct_singular(self, soln, iter, ierr, singular_index, reduced_from, reduced_to, times_singular) ! Try to correct the singular Jacobian matrix ! Arguments @@ -1876,6 +1886,7 @@ subroutine EqSolver_correct_singular(self, soln, iter, ierr, singular_index, red integer, intent(out), optional :: singular_index integer, intent(out), optional :: reduced_from integer, intent(out), optional :: reduced_to + integer, intent(in), optional :: times_singular ! Locals integer :: i, j, k ! Iterators @@ -2013,7 +2024,8 @@ subroutine EqSolver_correct_singular(self, soln, iter, ierr, singular_index, red end if ! Component-reduction fallback for persistent element-row singularities. - if (.not. made_change .and. ierr >= 1 .and. ierr <= ne .and. iter < 1 .and. & + if (.not. made_change .and. ierr >= 1 .and. ierr <= ne .and. & + (iter < 1 .or. (present(times_singular) .and. times_singular >= 4)) .and. & ne > 1 .and. .not. (self%ions .and. self%active_ions)) then call log_info("Reducing active element equations after singular restart on "// & trim(self%products%element_names(ierr))) @@ -2458,10 +2470,21 @@ subroutine EqSolver_solve(self, soln, type, state1, state2, reactant_weights, pa times_singular = times_singular + 1 if (times_singular > 8) then soln%converged = .false. - if (EqSolution_has_transport_seed(soln)) then - call EqSolution_restore_seed(soln) - else if (soln%T > 0.0d0 .and. soln%n > 0.0d0) then - call EqSolution_restore_transport_seed(soln) + if (num_reduced > 0 .and. soln%T > 0.0d0 .and. soln%n > 0.0d0) then + ! Retain the reduced-component state at the current + ! point after repeated singular matrices instead of + ! restoring the previous seed. + if (self%transport) then + call self%update_transport_basis(soln) + call compute_transport_properties(self, soln) + end if + call self%post_process(soln, .false.) + else + if (EqSolution_has_transport_seed(soln)) then + call EqSolution_restore_seed(soln) + else if (soln%T > 0.0d0 .and. soln%n > 0.0d0) then + call EqSolution_restore_transport_seed(soln) + end if end if call EqSolver_restore_reduced_elements(self, soln, num_reduced, reduced_from, reduced_to) call log_warning('EqSolver_solve: Too many singular matrices encountered.') @@ -2472,7 +2495,8 @@ subroutine EqSolver_solve(self, soln, type, state1, state2, reactant_weights, pa singular_index_iter = 0 reduced_from_iter = 0 reduced_to_iter = 0 - call self%correct_singular(soln, iter, ierr, singular_index_iter, reduced_from_iter, reduced_to_iter) + call self%correct_singular(soln, iter, ierr, singular_index_iter, reduced_from_iter, reduced_to_iter, & + times_singular) if (singular_index_iter > 0) singular_index = singular_index_iter if (reduced_to_iter > 0) then num_reduced = num_reduced + 1 @@ -5063,7 +5087,7 @@ subroutine compute_transport_properties(eq_solver, eq_soln, frozen_shock) ! Compute the transport properties of a mixture for a given equilibrium solution ! Arguments - class(EqSolver), target :: eq_solver + class(EqSolver), intent(in), target :: eq_solver type(EqSolution), intent(inout) :: eq_soln logical, intent(in), optional :: frozen_shock ! TODO: Estimate mole fractions for a frozen shock problem diff --git a/source/main.f90 b/source/main.f90 index 388cb56..851bdd8 100644 --- a/source/main.f90 +++ b/source/main.f90 @@ -6,7 +6,8 @@ program cea use cea_transport, only: TransportDB, read_transport use cea_input, only: ProblemDB, read_input use cea_equilibrium, only: EqSolver, EqSolution, EqPartials - use cea_rocket, only: RocketSolver, RocketSolution + use cea_rocket, only: RocketSolver, RocketSolution, & + rocket_status_partial, rocket_warning_condensed_temp_range use cea_shock, only: ShockSolver, ShockSolution use cea_detonation, only: DetonSolver, DetonSolution use cea_db_compile, only: compile_thermo_database, compile_transport_database @@ -400,7 +401,7 @@ subroutine run_rocket_problem(prob, thermo, solver, solutions) real(dp), allocatable :: weights(:) integer :: i, j, k, num_pc, num_of real(dp) :: pc, hc - real(dp), allocatable :: tc, tc_est, mdot, ac_at + real(dp) :: tc, tc_est, mdot, ac_at real(dp), allocatable :: subar(:), supar(:), pi_p(:) logical :: fac, frz, eql, need_hc integer :: nfrz @@ -408,6 +409,12 @@ subroutine run_rocket_problem(prob, thermo, solver, solutions) ! Initialize need_hc = .false. + hc = empty_dp + tc = empty_dp + tc_est = empty_dp + mdot = empty_dp + ac_at = empty_dp + allocate(pi_p(0), subar(0), supar(0)) ! Get the reactants Mixture object reactants = Mixture(thermo, input_reactants=prob%reactants, ions=prob%problem%include_ions) @@ -439,7 +446,7 @@ subroutine run_rocket_problem(prob, thermo, solver, solutions) end if ! Get the rocket variables - pi_p = prob%problem%pcp_schedule%values + if (allocated(prob%problem%pcp_schedule)) pi_p = prob%problem%pcp_schedule%values if (allocated(prob%problem%subar_schedule)) subar = prob%problem%subar_schedule%values if (allocated(prob%problem%supar_schedule)) supar = prob%problem%supar_schedule%values if (allocated(prob%problem%mdot)) mdot = prob%problem%mdot @@ -1907,7 +1914,8 @@ subroutine thermo_output(ioout, prob, solver, solutions, partials) do k = 1, size(solutions, 3) ! Loop over o/f ratio - ! Legacy-style section headers to align major output structure with CEA2. + ! Use the established section headers so the major output + ! structure stays aligned with the rest of the report. select case(prob%problem%type) case ("hp") section_title = "THERMODYNAMIC EQUILIBRIUM COMBUSTION PROPERTIES AT ASSIGNED" @@ -2253,6 +2261,9 @@ subroutine rocket_output(ioout, prob, solver, solutions) ! Get the number of stations ("np"), and the number of exit stations ("ne") np = solutions(i,j,k)%num_pts + if (solutions(i,j,k)%last_completed_idx > 0) then + np = min(np, solutions(i,j,k)%last_completed_idx) + end if ! Remove the solution at "inifnity" if this is an FAC problem if (prob%problem%rkt_finite_area) then @@ -2418,21 +2429,21 @@ subroutine rocket_output(ioout, prob, solver, solutions) write(ioout, thermo_fmt) " Visc, Millipoise", (solutions(i, j, k)%eq_soln(idx)%viscosity, idx=1,np) write(ioout, *) "" - ! Equilibrium properies - write(ioout, '(A)') " WITH EQUILIBRIUM REACTIONS" - if (prob%output%siunit) then - write(ioout, thermo_fmt) " Cp, kJ/(kg-K) ", (solutions(i, j, k)%eq_soln(idx)%cp_eq, idx=1,np) - write(ioout, thermo_fmt) " Conductivity ", & - (solutions(i, j, k)%eq_soln(idx)%conductivity_eq, idx=1,np) - else - write(ioout, thermo_fmt) " Cp, cal/(g-K) ", (solutions(i, j, k)%eq_soln(idx)%cp_eq/4.184d0, idx=1,np) - write(ioout, thermo_fmt) " Conductivity ", & - (solutions(i, j, k)%eq_soln(idx)%conductivity_eq/4.184d0, idx=1,np) + if (.not. frozen) then + write(ioout, '(A)') " WITH EQUILIBRIUM REACTIONS" + if (prob%output%siunit) then + write(ioout, thermo_fmt) " Cp, kJ/(kg-K) ", (solutions(i, j, k)%eq_soln(idx)%cp_eq, idx=1,np) + write(ioout, thermo_fmt) " Conductivity ", & + (solutions(i, j, k)%eq_soln(idx)%conductivity_eq, idx=1,np) + else + write(ioout, thermo_fmt) " Cp, cal/(g-K) ", (solutions(i, j, k)%eq_soln(idx)%cp_eq/4.184d0, idx=1,np) + write(ioout, thermo_fmt) " Conductivity ", & + (solutions(i, j, k)%eq_soln(idx)%conductivity_eq/4.184d0, idx=1,np) + end if + write(ioout, thermo_fmt) " Prandtl Number ", (solutions(i, j, k)%eq_soln(idx)%Pr_eq, idx=1,np) + write(ioout, *) "" end if - write(ioout, thermo_fmt) " Prandtl Number ", (solutions(i, j, k)%eq_soln(idx)%Pr_eq, idx=1,np) - write(ioout, *) "" - ! Frozen properties write(ioout, '(A)') " WITH FROZEN REACTIONS" if (prob%output%siunit) then write(ioout, thermo_fmt) " Cp, kJ/(kg-K) ", (solutions(i, j, k)%eq_soln(idx)%cp_fr, idx=1,np) @@ -2449,45 +2460,47 @@ subroutine rocket_output(ioout, prob, solver, solutions) end if ! Print the performance parameters - write(ioout, *) "" - write(ioout, '(A)') " PERFORMANCE PARAMETERS" - write(ioout, *) "" + if (np > 1) then + write(ioout, *) "" + write(ioout, '(A)') " PERFORMANCE PARAMETERS" + write(ioout, *) "" - perf_label = 'Ae/At ' - write(perf_fmt, '("(",i0,"(F13.4))")') np-1 - write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' - write(ioout, trim(perf_fmt)) (solutions(i, j, k)%ae_at(idx), idx=2,np) - if (prob%output%siunit) then - perf_label = 'C*, m/s ' - write(perf_fmt, '("(",i0,"(F13.2))")') np-1 - write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' - write(ioout, trim(perf_fmt)) (solutions(i, j, k)%c_star(idx), idx=2,np) - else - perf_label = 'C*, ft/s ' - write(perf_fmt, '("(",i0,"(F13.2))")') np-1 - write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' - write(ioout, trim(perf_fmt)) (solutions(i, j, k)%c_star(idx)*3.2808349d0, idx=2,np) - end if - perf_label = 'Cf ' - write(perf_fmt, '("(",i0,"(F13.4))")') np-1 - write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' - write(ioout, trim(perf_fmt)) (solutions(i, j, k)%cf(idx), idx=2,np) - if (prob%output%siunit) then - perf_label = 'Ivac, m/s ' - write(perf_fmt, '("(",i0,"(F13.2))")') np-1 - write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' - write(ioout, trim(perf_fmt)) (solutions(i, j, k)%i_vac(idx), idx=2,np) - perf_label = 'Isp, m/s ' - write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' - write(ioout, trim(perf_fmt)) (solutions(i, j, k)%i_sp(idx), idx=2,np) - else - perf_label = 'Ivac, lb-s/lb ' - write(perf_fmt, '("(",i0,"(F13.2))")') np-1 + perf_label = 'Ae/At ' + write(perf_fmt, '("(",i0,"(F13.4))")') np-1 write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' - write(ioout, trim(perf_fmt)) (solutions(i, j, k)%i_vac(idx)/9.80665d0, idx=2,np) - perf_label = 'Isp, lb-s/lb ' + write(ioout, trim(perf_fmt)) (solutions(i, j, k)%ae_at(idx), idx=2,np) + if (prob%output%siunit) then + perf_label = 'C*, m/s ' + write(perf_fmt, '("(",i0,"(F13.2))")') np-1 + write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' + write(ioout, trim(perf_fmt)) (solutions(i, j, k)%c_star(idx), idx=2,np) + else + perf_label = 'C*, ft/s ' + write(perf_fmt, '("(",i0,"(F13.2))")') np-1 + write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' + write(ioout, trim(perf_fmt)) (solutions(i, j, k)%c_star(idx)*3.2808349d0, idx=2,np) + end if + perf_label = 'Cf ' + write(perf_fmt, '("(",i0,"(F13.4))")') np-1 write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' - write(ioout, trim(perf_fmt)) (solutions(i, j, k)%i_sp(idx)/9.80665d0, idx=2,np) + write(ioout, trim(perf_fmt)) (solutions(i, j, k)%cf(idx), idx=2,np) + if (prob%output%siunit) then + perf_label = 'Ivac, m/s ' + write(perf_fmt, '("(",i0,"(F13.2))")') np-1 + write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' + write(ioout, trim(perf_fmt)) (solutions(i, j, k)%i_vac(idx), idx=2,np) + perf_label = 'Isp, m/s ' + write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' + write(ioout, trim(perf_fmt)) (solutions(i, j, k)%i_sp(idx), idx=2,np) + else + perf_label = 'Ivac, lb-s/lb ' + write(perf_fmt, '("(",i0,"(F13.2))")') np-1 + write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' + write(ioout, trim(perf_fmt)) (solutions(i, j, k)%i_vac(idx)/9.80665d0, idx=2,np) + perf_label = 'Isp, lb-s/lb ' + write(ioout, '(1x,A,1x,A13)', advance='no') perf_label, '' + write(ioout, trim(perf_fmt)) (solutions(i, j, k)%i_sp(idx)/9.80665d0, idx=2,np) + end if end if ! Set the trace output value @@ -2555,15 +2568,18 @@ subroutine rocket_output(ioout, prob, solver, solutions) write(ioout, '(A)') "PRODUCTS WHICH WERE CONSIDERED BUT WHOSE "// mass_or_mole //" FRACTIONS" write(ioout, '(A, 1PE13.6 ,A)') "WERE LESS THAN", trace, " FOR ALL ASSIGNED CONDITIONS" write(ioout, '(A)') "" - ncols = 5 - nrows = (num_trace-1)/ncols + 1 - do ii = 1, nrows - last_row_cols = merge(mod(num_trace, ncols), ncols, mod(num_trace, ncols) /= 0 .and. ii == nrows) - write(ioout, '(5A16)') (trace_names((ii-1)*ncols + jj), jj=1,last_row_cols) - end do + if (num_trace > 0) then + ncols = 5 + nrows = (num_trace-1)/ncols + 1 + do ii = 1, nrows + last_row_cols = merge(mod(num_trace, ncols), ncols, mod(num_trace, ncols) /= 0 .and. ii == nrows) + write(ioout, '(5A16)') (trace_names((ii-1)*ncols + jj), jj=1,last_row_cols) + end do + end if write(ioout, '(A)') "" write(ioout, '(A)') "NOTE. WEIGHT FRACTION OF FUEL IN TOTAL FUELS AND OF OXIDANT IN TOTAL OXIDANTS" write(ioout, '(A)') "" + call write_rocket_legacy_warning(ioout, solutions(i,j,k)) ! Write out the additional exit conditions if any are in overflow if (exit_extra > 0) then @@ -2699,31 +2715,31 @@ subroutine rocket_output(ioout, prob, solver, solutions) write(ioout, '(F13.3)') (solutions(i, j, k)%eq_soln(idx)%viscosity, idx=x,y) write(ioout, *) "" - ! Equilibrium properies - write(ioout, '(A)') " WITH EQUILIBRIUM REACTIONS" - if (prob%output%siunit) then - write(ioout, '(A, F13.4)', advance="no") " Cp, kJ/(kg-K) ", & - (solutions(i, j, k)%eq_soln(idx)%cp_eq, idx=1,nc) - write(ioout, '(F13.4)') (solutions(i, j, k)%eq_soln(idx)%cp_eq, idx=x,y) + if (.not. frozen) then + write(ioout, '(A)') " WITH EQUILIBRIUM REACTIONS" + if (prob%output%siunit) then + write(ioout, '(A, F13.4)', advance="no") " Cp, kJ/(kg-K) ", & + (solutions(i, j, k)%eq_soln(idx)%cp_eq, idx=1,nc) + write(ioout, '(F13.4)') (solutions(i, j, k)%eq_soln(idx)%cp_eq, idx=x,y) - write(ioout, '(A, F13.3)', advance="no") " Conductivity ", & - (solutions(i, j, k)%eq_soln(idx)%conductivity_eq, idx=1,nc) - write(ioout, '(F13.3)') (solutions(i, j, k)%eq_soln(idx)%conductivity_eq, idx=x,y) - else - write(ioout, '(A, F13.4)', advance="no") " Cp, cal/(g-K) ", & - (solutions(i, j, k)%eq_soln(idx)%cp_eq/4.184d0, idx=1,nc) - write(ioout, '(F13.4)') (solutions(i, j, k)%eq_soln(idx)%cp_eq/4.184d0, idx=x,y) + write(ioout, '(A, F13.3)', advance="no") " Conductivity ", & + (solutions(i, j, k)%eq_soln(idx)%conductivity_eq, idx=1,nc) + write(ioout, '(F13.3)') (solutions(i, j, k)%eq_soln(idx)%conductivity_eq, idx=x,y) + else + write(ioout, '(A, F13.4)', advance="no") " Cp, cal/(g-K) ", & + (solutions(i, j, k)%eq_soln(idx)%cp_eq/4.184d0, idx=1,nc) + write(ioout, '(F13.4)') (solutions(i, j, k)%eq_soln(idx)%cp_eq/4.184d0, idx=x,y) - write(ioout, '(A, F13.3)', advance="no") " Conductivity ", & - (solutions(i, j, k)%eq_soln(idx)%conductivity_eq/4.184d0, idx=1,nc) - write(ioout, '(F13.3)') (solutions(i, j, k)%eq_soln(idx)%conductivity_eq/4.184d0, idx=x,y) + write(ioout, '(A, F13.3)', advance="no") " Conductivity ", & + (solutions(i, j, k)%eq_soln(idx)%conductivity_eq/4.184d0, idx=1,nc) + write(ioout, '(F13.3)') (solutions(i, j, k)%eq_soln(idx)%conductivity_eq/4.184d0, idx=x,y) + end if + write(ioout, '(A, F13.4)', advance="no") " Prandtl Number ", & + (solutions(i, j, k)%eq_soln(idx)%Pr_eq, idx=1,nc) + write(ioout, '(F13.4)') (solutions(i, j, k)%eq_soln(idx)%Pr_eq, idx=x,y) + write(ioout, *) "" end if - write(ioout, '(A, F13.4)', advance="no") " Prandtl Number ", & - (solutions(i, j, k)%eq_soln(idx)%Pr_eq, idx=1,nc) - write(ioout, '(F13.4)') (solutions(i, j, k)%eq_soln(idx)%Pr_eq, idx=x,y) - write(ioout, *) "" - ! Frozen properties write(ioout, '(A)') " WITH FROZEN REACTIONS" if (prob%output%siunit) then write(ioout, '(A, F13.4)', advance="no") " Cp, kJ/(kg-K) ", & @@ -2866,6 +2882,20 @@ subroutine rocket_output(ioout, prob, solver, solutions) end subroutine + subroutine write_rocket_legacy_warning(ioout, solution) + integer, intent(in) :: ioout + type(RocketSolution), intent(in) :: solution + + if (solution%status_code /= rocket_status_partial) return + + select case (solution%warning_code) + case (rocket_warning_condensed_temp_range) + write(ioout, '(A)') " WARNING!! CALCULATIONS WERE STOPPED BECAUSE NEXT POINT IS MORE" + write(ioout, '(A)') " THAN 50 K BELOW THE TEMPERATURE RANGE OF A CONDENSED SPECIES (ROCKET)" + write(ioout, '(A)') "" + end select + end subroutine + subroutine compute_fuel_ratios(prob, reactants, idx, of_ratio, pct_fuel, r_eq, phi_eq) ! Compute all of the fuel ratio values (o/f ratio, % fuel, r eq. ratio, phi eq. ratio) at once @@ -2883,9 +2913,13 @@ subroutine compute_fuel_ratios(prob, reactants, idx, of_ratio, pct_fuel, r_eq, p real(dp) :: ratio_val real(dp), allocatable :: fuel_moles(:), oxidant_moles(:) real(dp), allocatable :: fuel_weights(:), oxidant_weights(:) + real(dp), allocatable :: vm(:), vp(:), b0_ox(:), b0_fu(:) + real(dp) :: vm_fu, vp_fu, vm_ox, vp_ox allocate(fuel_moles(reactants%num_species), oxidant_moles(reactants%num_species), & - fuel_weights(reactants%num_species), oxidant_weights(reactants%num_species)) + fuel_weights(reactants%num_species), oxidant_weights(reactants%num_species), & + vm(reactants%num_elements), vp(reactants%num_elements), & + b0_ox(reactants%num_elements), b0_fu(reactants%num_elements)) fuel_weights = 0.0d0 oxidant_weights = 0.0d0 @@ -2993,6 +3027,24 @@ subroutine compute_fuel_ratios(prob, reactants, idx, of_ratio, pct_fuel, r_eq, p end if + if (sum(fuel_weights) > 0.0d0 .and. sum(oxidant_weights) > 0.0d0 .and. of_ratio > 0.0d0) then + call reactants%get_valence(vm, vp) + b0_ox = reactants%element_amounts_from_weights(oxidant_weights) + b0_fu = reactants%element_amounts_from_weights(fuel_weights) + + vm_ox = dot_product(vm, b0_ox) + vp_ox = dot_product(vp, b0_ox) + vm_fu = dot_product(vm, b0_fu) + vp_fu = dot_product(vp, b0_fu) + + if (abs(vm_ox + vp_ox) > tiny(1.0d0)) then + phi_eq = reactants%phi_from_of(oxidant_weights, fuel_weights, of_ratio) + end if + if (abs(vm_fu + vm_ox*of_ratio) > tiny(1.0d0)) then + r_eq = reactants%equivalence_from_of(oxidant_weights, fuel_weights, of_ratio) + end if + end if + end if end subroutine diff --git a/source/rocket.f90 b/source/rocket.f90 index 1233c25..b127f04 100644 --- a/source/rocket.f90 +++ b/source/rocket.f90 @@ -10,6 +10,13 @@ module cea_rocket use fb_utils implicit none + integer, parameter :: rocket_status_success = 0 + integer, parameter :: rocket_status_partial = 1 + integer, parameter :: rocket_status_failed = 2 + + integer, parameter :: rocket_warning_none = 0 + integer, parameter :: rocket_warning_condensed_temp_range = 1 + type :: RocketSolver !! Rocket solver class @@ -85,6 +92,12 @@ module cea_rocket ! Convergence variables logical :: converged = .false. !! Convergence flag + integer :: status_code = rocket_status_failed + !! Rocket solve status + integer :: warning_code = rocket_warning_none + !! Legacy warning selector for partial output + integer :: last_completed_idx = 0 + !! Last successfully completed station index end type interface RocketSolution @@ -93,6 +106,42 @@ module cea_rocket contains + logical function has_real_value(value) + real(dp), intent(in), optional :: value + + has_real_value = .false. + if (.not. present(value)) return + has_real_value = (value /= empty_dp) + end function + + logical function has_array_values(values) + real(dp), intent(in), optional :: values(:) + + has_array_values = .false. + if (.not. present(values)) return + has_array_values = (size(values) > 0) + end function + + subroutine frozen_fit_bounds(species, T_ref, T_low, T_high) + use cea_thermo, only: SpeciesThermo + type(SpeciesThermo), intent(in) :: species + real(dp), intent(in) :: T_ref + real(dp), intent(out) :: T_low + real(dp), intent(out) :: T_high + + integer :: i, idx + + idx = 1 + do i = 1, species%num_intervals + if (T_ref > species%T_fit(i, 1)) then + idx = i + end if + end do + + T_low = species%T_fit(idx, 1) + T_high = species%T_fit(idx, 2) + end subroutine + !----------------------------------------------------------------------- ! RocketSolver !----------------------------------------------------------------------- @@ -140,7 +189,8 @@ subroutine RocketSolver_frozen(self, soln, idx, n_frz) real(dp) :: cpsum, ssum ! Temporary variables for mixture properties real(dp) :: cpj, sj ! Temporary variables for species properties real(dp) :: dlnt ! Update variable for log-temperature - real(dp) :: T_low, T_high ! Condensed species temperature bounds [K] + real(dp) :: T_low, T_high ! Species temperature bounds [K] + real(dp) :: gas_T_min ! Minimum valid gas-fit lower bound [K] call log_debug("Starting frozen calculations") @@ -189,15 +239,27 @@ subroutine RocketSolver_frozen(self, soln, idx, n_frz) soln%eq_partials(idx)%dlnV_dlnT = 1.0d0 call self%eq_solver%products%calc_thermo(soln%eq_soln(idx)%thermo, soln%eq_soln(idx)%T) - ! Stop if any frozen condensed species is outside its - ! valid temperature range by more than 50 K. + gas_T_min = huge(1.0d0) + do j = 1, ng + if (abs(soln%eq_soln(n_frz)%nj(j)) <= approx_zero_tol) cycle + gas_T_min = min(gas_T_min, minval(self%eq_solver%products%species(j)%T_fit(:, 1))) + end do + + if (soln%eq_soln(idx)%T < (0.8d0 * gas_T_min)) then + call log_warning("Frozen calculations stopped: temperature below lower-limit "// & + "(0.8*Tmin="//to_str(0.8d0 * gas_T_min)//")") + call mark_partial_stop(soln, idx-1, rocket_warning_condensed_temp_range) + return + end if + + ! Check the thermo interval associated with the frozen + ! composition point, not the species' overall span. in_range = .true. do j = 1, self%eq_solver%num_condensed ! TODO(smooth_truncation): smooth gating means species are rarely exactly zero. ! Frozen-mode checks intentionally use a practical-zero tolerance. if (abs(soln%eq_soln(n_frz)%nj(ng+j)) <= approx_zero_tol) cycle - T_low = minval(self%eq_solver%products%species(ng+j)%T_fit(:, 1)) - T_high = maxval(self%eq_solver%products%species(ng+j)%T_fit(:, 2)) + call frozen_fit_bounds(self%eq_solver%products%species(ng+j), soln%eq_soln(n_frz)%T, T_low, T_high) if (soln%eq_soln(idx)%T < (T_low-phase_gap) .or. soln%eq_soln(idx)%T > (T_high+phase_gap)) then in_range = .false. exit @@ -207,7 +269,7 @@ subroutine RocketSolver_frozen(self, soln, idx, n_frz) if (.not. in_range) then call log_warning("Frozen calculations stopped: temperature is more than 50 K outside "// & "the range of a condensed species") - soln%converged = .false. + call mark_partial_stop(soln, idx-1, rocket_warning_condensed_temp_range) return end if @@ -223,7 +285,12 @@ subroutine RocketSolver_frozen(self, soln, idx, n_frz) if (.not. finalized) then call log_warning("Frozen calculations did not converge in 8 iterations") - soln%converged = .false. + if (idx > 1) then + call mark_partial_stop(soln, idx-1, rocket_warning_condensed_temp_range) + else + soln%converged = .false. + soln%status_code = rocket_status_failed + end if return end if @@ -237,7 +304,28 @@ subroutine RocketSolver_frozen(self, soln, idx, n_frz) soln%eq_soln(idx)%energy = soln%eq_soln(idx)%enthalpy - soln%eq_soln(idx)%n*soln%eq_soln(idx)%T*R/1.d3 soln%eq_soln(idx)%entropy = soln%eq_soln(n_frz)%entropy soln%eq_soln(idx)%gibbs_energy = (soln%eq_soln(idx)%enthalpy - soln%eq_soln(idx)%T*soln%eq_soln(idx)%entropy) + call self%eq_solver%compute_transport_state(soln%eq_soln(idx)) + end subroutine + subroutine mark_completed(soln, idx) + type(RocketSolution), intent(inout) :: soln + integer, intent(in) :: idx + + soln%last_completed_idx = max(soln%last_completed_idx, idx) + soln%converged = .true. + soln%status_code = rocket_status_success + end subroutine + + subroutine mark_partial_stop(soln, last_idx, warning_code) + type(RocketSolution), intent(inout) :: soln + integer, intent(in) :: last_idx + integer, intent(in) :: warning_code + + soln%converged = .false. + soln%status_code = rocket_status_partial + soln%warning_code = warning_code + soln%last_completed_idx = max(soln%last_completed_idx, last_idx) + soln%num_pts = max(0, soln%last_completed_idx) end subroutine subroutine RocketSolver_solve_throat(self, soln, idx, pc, h_inf, s0, weights, awt) @@ -393,6 +481,8 @@ subroutine RocketSolver_solve_throat_frozen(self, soln, idx, n_frz, pc, h_inf, a end do + call mark_completed(soln, idx) + end subroutine subroutine RocketSolver_solve_pi_p(self, soln, idx, pc, pi_p, h_inf, s0, weights) @@ -453,6 +543,7 @@ subroutine RocketSolver_solve_pi_p(self, soln, idx, pc, pi_p, h_inf, s0, weights ! Iterate idx = idx + 1 + call mark_completed(soln, idx-1) end do @@ -484,8 +575,8 @@ subroutine RocketSolver_solve_pi_p_frozen(self, soln, idx, n_frz, pc, pi_p, h_in (soln%pressure(soln%throat_idx)*soln%v_sonic(soln%throat_idx)) pip_nf = pc/soln%pressure(n_frz) do i = 1, size(pi_p) - ! Legacy frozen scheduling: omit assigned pressure ratios lower than - ! the value at the freeze point. + ! Omit assigned pressure ratios lower than the value at the + ! freeze point. if (pi_p(i) < pip_nf) then call log_info('RocketSolver: WARNING!! FOR FROZEN PERFORMANCE, POINT OMITTED BECAUSE '// & 'ASSIGNED pi/p IS LESS THAN VALUE AT nfz='//to_str(n_frz)) @@ -514,6 +605,7 @@ subroutine RocketSolver_solve_pi_p_frozen(self, soln, idx, n_frz, pc, pi_p, h_in soln%ae_at(idx) = soln%eq_soln(idx)%n*soln%eq_soln(idx)%T/(soln%pressure(idx)*usq**0.5*awt) idx = idx + 1 + call mark_completed(soln, idx-1) end do @@ -603,6 +695,7 @@ subroutine RocketSolver_solve_subar(self, soln, idx, pc, subar, h_inf, s0, weigh ! Iterate to the next exit condition idx = idx + 1 + call mark_completed(soln, idx-1) end do @@ -688,6 +781,7 @@ subroutine RocketSolver_solve_supar(self, soln, idx, pc, supar, h_inf, s0, weigh end if idx = idx + 1 + call mark_completed(soln, idx-1) end do @@ -776,6 +870,7 @@ subroutine RocketSolver_solve_supar_frozen(self, soln, idx, n_frz, pc, supar, h_ end do idx = idx + 1 + call mark_completed(soln, idx-1) end do @@ -822,9 +917,10 @@ subroutine RocketSolver_solve_iac(self, soln, reactant_weights, pc, pi_p, subar, call log_debug("Initialized RocketSolution") soln%throat_idx = 2 soln%converged = .true. + soln%status_code = rocket_status_success ! Set the equilibrium problem type - if (present(hc)) then + if (has_real_value(hc)) then prob_type = "hp" state1 = hc h_inf = hc @@ -851,7 +947,7 @@ subroutine RocketSolver_solve_iac(self, soln, reactant_weights, pc, pi_p, subar, ! ----------------------------------------------- ! Chamber conditions (infinity) ! ----------------------------------------------- - if (present(tc_est)) then + if (has_real_value(tc_est)) then soln%eq_soln(1) = EqSolution(self%eq_solver, T_init=tc_est) else soln%eq_soln(1) = EqSolution(self%eq_solver) @@ -879,6 +975,7 @@ subroutine RocketSolver_solve_iac(self, soln, reactant_weights, pc, pi_p, subar, soln%eq_soln(1)%gamma_s = soln%eq_partials(1)%gamma_s soln%v_sonic(1) = sqrt(soln%eq_soln(1)%n*R*soln%gamma_s(1)*soln%eq_soln(1)%T) end if + call mark_completed(soln, 1) ! Save some chamber solution variables for later use state1 = soln%eq_soln(1)%calc_entropy_sum(self%eq_solver) ! Combustor entropy @@ -894,13 +991,17 @@ subroutine RocketSolver_solve_iac(self, soln, reactant_weights, pc, pi_p, subar, else call self%solve_throat(soln, idx, pc, h_inf, state1, reactant_weights, awt) end if - if (.not. soln%converged) return + if (soln%status_code == rocket_status_failed) return + if (soln%status_code == rocket_status_partial) then + call self%post_process(soln, .false.) + return + end if ln_pinf_pt = log(soln%pressure(1)/soln%pressure(2)) ! ----------------------------------------------- ! Exit conditions: pressure ratio ! ----------------------------------------------- - if (present(pi_p)) then + if (has_array_values(pi_p)) then idx = 3 if (frozen .and. idx > n_frz_) then @@ -908,7 +1009,11 @@ subroutine RocketSolver_solve_iac(self, soln, reactant_weights, pc, pi_p, subar, else call self%solve_pi_p(soln, idx, pc, pi_p, h_inf, state1, reactant_weights) end if - if (.not. soln%converged) return + if (soln%status_code == rocket_status_failed) return + if (soln%status_code == rocket_status_partial) then + call self%post_process(soln, .false.) + return + end if else ! If pi_p not present, idx stays at 3 for subsequent sections idx = 3 @@ -918,7 +1023,7 @@ subroutine RocketSolver_solve_iac(self, soln, reactant_weights, pc, pi_p, subar, ! Exit conditions: subsonic area ratio ! ----------------------------------------------- - if (present(subar)) then + if (has_array_values(subar)) then if (frozen .and. n_frz_ > 1) then call log_info('RocketSolver: WARNING!! FREEZING IS NOT ALLOWED AT A SUBSONIC PRESSURE RATIO') @@ -932,14 +1037,18 @@ subroutine RocketSolver_solve_iac(self, soln, reactant_weights, pc, pi_p, subar, ! Exit conditions: supersonic area ratio ! ----------------------------------------------- - if (present(supar)) then + if (has_array_values(supar)) then if (frozen .and. idx > n_frz_) then call self%solve_supar_frozen(soln, idx, n_frz, pc, supar, h_inf, reactant_weights, idx-1, ln_pinf_pt, awt) else call self%solve_supar(soln, idx, pc, supar, h_inf, state1, reactant_weights, ln_pinf_pt, awt) end if - if (.not. soln%converged) return + if (soln%status_code == rocket_status_failed) return + if (soln%status_code == rocket_status_partial) then + call self%post_process(soln, .false.) + return + end if end if @@ -1009,28 +1118,29 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar, ! Set the total number of evaluation points num_pts = 4 ! injector + infinity + combustor + throat - if (present(pi_p)) num_pts = num_pts + size(pi_p) - if (present(subar)) num_pts = num_pts + size(subar) - if (present(supar)) num_pts = num_pts + size(supar) + if (has_array_values(pi_p)) num_pts = num_pts + size(pi_p) + if (has_array_values(subar)) num_pts = num_pts + size(subar) + if (has_array_values(supar)) num_pts = num_pts + size(supar) ! Initialize the solution variables soln = RocketSolution(self, num_pts=num_pts) call log_debug("Initialized RocketSolution") soln%throat_idx = 4 soln%converged = .true. + soln%status_code = rocket_status_success ! Set a flag to determine whether to use the contraction ratio or mdot for the chamber solution - if (present(ac_at)) then + if (has_real_value(ac_at)) then use_acat = .true. ac_at_ = ac_at - else if (present(mdot)) then + else if (has_real_value(mdot)) then use_acat = .false. else call abort("Either Ac/At or mdot/At must be specified for FAC rocket problem") end if ! Set the equilibrium problem type - if (present(hc)) then + if (has_real_value(hc)) then prob_type = "hp" state1 = hc h_inj = hc @@ -1060,7 +1170,7 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar, idx = 1 soln%station(idx) = "injector" - if (present(tc_est)) then + if (has_real_value(tc_est)) then soln%eq_soln(idx) = EqSolution(self%eq_solver, T_init=tc_est) soln%eq_soln(2) = EqSolution(self%eq_solver, T_init=tc_est) ! Initialize soln at infinity too else @@ -1076,6 +1186,7 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar, gamma_s = soln%eq_partials(idx)%gamma_s ! Save gamma_s as the initial value for the next station soln%gamma_s(1) = gamma_s soln%v_sonic(idx) = sqrt(soln%eq_soln(idx)%n*R*soln%gamma_s(idx)*soln%eq_soln(idx)%T) + call mark_completed(soln, idx) ! Save the injector enthalpy to compute conditions at infinity, which is isenthropic with the injector if (prob_type == "tp") h_inj = dot_product(soln%eq_soln(idx)%nj, soln%eq_soln(idx)%thermo%enthalpy)*soln%eq_soln(idx)%T @@ -1095,7 +1206,8 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar, ! P_c initial guess soln%pressure(3) = p_inf - ! Iterate until combustor conditions converge (legacy CEA2-style stopping test). + ! Iterate until combustor conditions converge with the same stopping + ! test used for this chamber closure loop. chamber_iter = 0 do chamber_iter = chamber_iter + 1 @@ -1112,8 +1224,8 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar, ! Throat conditions ! ----------------------------------------------- idx = 4 - ! Legacy SETEN-style initialization for FAC throat: - ! seed the throat solve from point 2 (infinity) by saving/using it. + ! Seed the throat solve from point 2 (infinity) by saving and + ! reusing that state as the initial condition. soln%eq_soln(3) = soln%eq_soln(2) soln%i_save = -2 call self%solve_throat(soln, idx, p_inf, h_inj, S_ref, reactant_weights, awt) @@ -1224,12 +1336,16 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar, if (frozen .and. idx > n_frz_) then call self%solve_throat_frozen(soln, idx, n_frz_, p_inf, h_inj, awt) end if - if (.not. soln%converged) return + if (soln%status_code == rocket_status_failed) return + if (soln%status_code == rocket_status_partial) then + call self%post_process(soln, .true.) + return + end if ! ----------------------------------------------- ! Exit conditions: pressure ratio ! ----------------------------------------------- - if (present(pi_p)) then + if (has_array_values(pi_p)) then idx = 5 if (frozen .and. idx > n_frz_) then @@ -1237,7 +1353,11 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar, else call self%solve_pi_p(soln, idx, pc, pi_p, h_inj, S_ref, reactant_weights) end if - if (.not. soln%converged) return + if (soln%status_code == rocket_status_failed) return + if (soln%status_code == rocket_status_partial) then + call self%post_process(soln, .true.) + return + end if else ! If pi_p not present, idx stays at 5 for subsequent sections idx = 5 @@ -1247,7 +1367,7 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar, ! Exit conditions: subsonic area ratio ! ----------------------------------------------- - if (present(subar)) then + if (has_array_values(subar)) then if (frozen .and. n_frz_ > 1) then call log_info('RocketSolver: WARNING!! FREEZING IS NOT ALLOWED AT A SUBSONIC PRESSURE RATIO') @@ -1263,7 +1383,7 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar, ! Exit conditions: supersonic area ratio ! ----------------------------------------------- - if (present(supar)) then + if (has_array_values(supar)) then ln_pinf_pt = log(soln%pressure(2)/soln%pressure(4)) if (frozen .and. idx > n_frz_) then @@ -1273,7 +1393,11 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar, call self%solve_supar(soln, idx, soln%pressure(2), supar, h_inj, S_ref, reactant_weights, & ln_pinf_pt, awt) end if - if (.not. soln%converged) return + if (soln%status_code == rocket_status_failed) return + if (soln%status_code == rocket_status_partial) then + call self%post_process(soln, .true.) + return + end if end if ! Omitted frozen schedule points do not consume output indices. diff --git a/source/rocket_test.pf b/source/rocket_test.pf index 219fcac..c8b4a01 100644 --- a/source/rocket_test.pf +++ b/source/rocket_test.pf @@ -464,6 +464,83 @@ contains end subroutine + @test + subroutine test_rocket_frozen_partial_output_fail1() + type(Mixture) :: products + type(Mixture) :: reactants + type(RocketSolver) :: solver + type(RocketSolution) :: soln + character(:), allocatable :: product_names(:) + real(dp) :: hc, pc, weights(4) + real(dp), parameter :: tol = 1.0d-4 + + reactants = Mixture(all_thermo, ['CH4(L) ', 'C2H6(L)', 'C3H8(L)', 'O2(L) ']) + product_names = reactants%get_products(all_thermo) + products = Mixture(all_thermo, product_names) + + weights = reactants%weights_from_of([0.0d0, 0.0d0, 0.0d0, 1.0d0], & + [0.95d0, 0.036d0, 0.014d0, 0.0d0], 0.7777778d0) + hc = -150.411d0 + pc = psi_to_bar(1000.0d0) + + solver = RocketSolver(products, reactants) + soln = solver%solve(weights, pc, supar=[10.433d0], hc=hc, fac=.false., n_frz=1) + + @assertFalse(soln%converged) + @assertEqual(rocket_status_partial, soln%status_code) + @assertEqual(rocket_warning_condensed_temp_range, soln%warning_code) + @assertEqual(2, soln%num_pts) + @assertEqual(2, soln%last_completed_idx) + + @assertRelativelyEqual(1369.34d0, soln%eq_soln(1)%T, tol) + @assertRelativelyEqual(68.947d0, soln%pressure(1), tol) + @assertRelativelyEqual(1210.78d0, soln%eq_soln(2)%T, tol) + @assertRelativelyEqual(38.038d0, soln%pressure(2), tol) + @assertRelativelyEqual(1.0d0, soln%mach(2), tol) + @assertRelativelyEqual(1466.5d0, soln%c_star(2), tol) + end subroutine + + @test + subroutine test_rocket_frozen_partial_output_fail2() + type(Mixture) :: products + type(Mixture) :: reactants + type(RocketSolver) :: solver + type(RocketSolution) :: soln + real(dp) :: hc, pc, weights(2) + real(dp), parameter :: tol = 1.0d-4 + real(dp), parameter :: transport_tol = 5.0d-4 + + reactants = Mixture(all_thermo, ['CH4', 'N2 ']) + products = Mixture(all_thermo, ['CH4', 'N2 ']) + + weights = reactants%weights_from_of([0.0d0, 1.0d0], [1.0d0, 0.0d0], 13.286d0) + hc = reactants%calc_enthalpy(weights, [1000.0d0, 1000.0d0])/R + pc = psi_to_bar(800.0d0) + + solver = RocketSolver(products, reactants, all_transport=all_transport) + soln = solver%solve(weights, pc, & + supar=[10.0d0, 20.0d0, 30.0d0, 40.0d0, 50.0d0, 60.0d0, 70.0d0, 80.0d0, 90.0d0, 100.0d0, 120.0d0, 135.0d0, 150.0d0], & + hc=hc, fac=.false., n_frz=1) + + @assertFalse(soln%converged) + @assertEqual(rocket_status_partial, soln%status_code) + @assertEqual(rocket_warning_condensed_temp_range, soln%warning_code) + @assertEqual(6, soln%num_pts) + @assertEqual(6, soln%last_completed_idx) + + @assertRelativelyEqual(1000.0d0, soln%eq_soln(1)%T, tol) + @assertRelativelyEqual(871.91d0, soln%eq_soln(2)%T, tol) + @assertRelativelyEqual(294.16d0, soln%eq_soln(3)%T, tol) + @assertRelativelyEqual(220.95d0, soln%eq_soln(4)%T, tol) + @assertRelativelyEqual(187.22d0, soln%eq_soln(5)%T, tol) + @assertRelativelyEqual(166.58d0, soln%eq_soln(6)%T, tol) + @assertRelativelyEqual(40.0d0, soln%ae_at(6), tol) + @assertRelativelyEqual(839.1d0, soln%c_star(6), tol) + @assertRelativelyEqual(1.7089d0, soln%cf(6), tol) + @assertRelativelyEqual(1.1219d0, soln%eq_soln(3)%cp_fr, tol) + @assertRelativelyEqual(0.2654d0, soln%eq_soln(3)%conductivity_fr, transport_tol) + end subroutine + @test subroutine test_rocket_transport type(Mixture) :: products