initial fix for w5 velocity error

Author: markleader

source/equilibrium_test.pf

@@ -470,6 +470,17 @@ contains
         ! Test the full solve
         call solver%solve(solution, 'uv', u_reac, v_reac, weights)
 
+        ! write(*,*) "T = ", solution%T, " rel = ", &
+        !    abs((2418.5382574474097d0 - solution%T)/2418.5382574474097d0)
+        ! write(*,*) "ln_nj(1) = ", solution%ln_nj(1), " rel = ", &
+        !    abs((-52.374890133250204d0 - solution%ln_nj(1))/(-52.374890133250204d0))
+        ! write(*,*) "ln_nj(2) = ", solution%ln_nj(2), " rel = ", &
+        !    abs((-8.0940265640663025d0 - solution%ln_nj(2))/(-8.0940265640663025d0))
+        ! write(*,*) "ln_nj(3) = ", solution%ln_nj(3), " rel = ", &
+        !    abs((-42.957051147438705d0 - solution%ln_nj(3))/(-42.957051147438705d0))
+        ! write(*,*) "ln_nj(7) = ", solution%ln_nj(7), " rel = ", &
+        !    abs((-9.7591561204275976d0 - solution%ln_nj(7))/(-9.7591561204275976d0))
+
         @assertRelativelyEqual(2418.5382574474097d0, solution%T, tol)
         @assertRelativelyEqual(-52.374890133250204d0, solution%ln_nj(1), tol)
         @assertRelativelyEqual(-8.0940265640663025d0, solution%ln_nj(2), tol)

source/input_test.pf

@@ -9,7 +9,6 @@ contains
     subroutine test_read_input
         type(ProblemDB), allocatable :: problems(:)
         allocate(problems(0))
-        write(*,*) 'Testing read_input'
         !call set_log_level(log_levels%debug)
         problems = read_input('samples/rp1311_examples.inp')
         @assertEqual(14, size(problems))

source/main.f90

@@ -872,7 +872,7 @@ subroutine shock_output(ioout, prob, solver, solutions)
         real(dp) :: trace
         logical :: incident, reflected, equilibrium, frozen, input_reflected
         logical :: write_incd_frz, write_refl_frz, write_incd_eql, write_refl_eql
-        logical :: use_mach
+        logical :: use_mach, have_incident_state, have_reflected_state
         character(snl), allocatable :: trace_names(:)
         logical, allocatable :: is_trace(:)
         character(:), allocatable :: eq_fmt
@@ -945,7 +945,7 @@ subroutine shock_output(ioout, prob, solver, solutions)
                 incd_type = "FROZEN     "
                 refl_type = "FROZEN     "
                 write_refl_eql = .false.
-                if (reflected) then
+                if (input_reflected) then
                     write_refl_frz = .true.
                 else
                     write_refl_frz = .false.
@@ -1021,6 +1021,22 @@ subroutine shock_output(ioout, prob, solver, solutions)
                 end if
             end if
 
+            have_incident_state = .false.
+            have_reflected_state = .false.
+            do i = 1, m
+                if (solutions(i, 1, k)%eq_soln(2)%T > 0.0d0) have_incident_state = .true.
+                if (solutions(i, 1, k)%eq_soln(3)%T > 0.0d0) have_reflected_state = .true.
+            end do
+            if (.not. have_incident_state) then
+                write_incd_frz = .false.
+                write_incd_eql = .false.
+                write_refl_frz = .false.
+                write_refl_eql = .false.
+            else if (.not. have_reflected_state) then
+                write_refl_frz = .false.
+                write_refl_eql = .false.
+            end if
+
             ! -------------------------------------------------------------------
             ! Write results for unshocked gas (equilibrium)
             ! -------------------------------------------------------------------

source/rocket_test.pf

@@ -43,6 +43,10 @@ contains
         ! Products order:  H, H2, H2O, H2O2, HO2, O, O2, O3, OH, H2O(L), H2O(cr)
 
         ! Test the chamber conditions
+        ! write(*,*) "st1 T = ", soln%eq_soln(1)%T, " rel = ", &
+        !    abs((3383.8446259451043d0 - soln%eq_soln(1)%T)/3383.8446259451043d0)
+        ! write(*,*) "st1 ln_nj(H2) = ", soln%eq_soln(1)%ln_nj(2), " rel = ", &
+        !    abs((-3.7643107135196541d0 - soln%eq_soln(1)%ln_nj(2))/(-3.7643107135196541d0))
         @assertRelativelyEqual(3383.8446259451043d0,  soln%eq_soln(1)%T, tol)
         @assertRelativelyEqual(-5.9391066300902899d0, soln%eq_soln(1)%ln_nj(1), tol)
         @assertRelativelyEqual(-3.7643107135196541d0, soln%eq_soln(1)%ln_nj(2), tol)
@@ -58,6 +62,10 @@ contains
         @assertRelativelyEqual(1.1446952882630981d0, soln%eq_partials(1)%gamma_s, tol)
 
         ! Test the throat conditions
+        ! write(*,*) "st2 T = ", soln%eq_soln(2)%T, " rel = ", &
+        !    abs((3185.6731730807096d0 - soln%eq_soln(2)%T)/3185.6731730807096d0)
+        ! write(*,*) "st2 ln_nj(H2) = ", soln%eq_soln(2)%ln_nj(2), " rel = ", &
+        !    abs((-3.7758918668469459d0 - soln%eq_soln(2)%ln_nj(2))/(-3.7758918668469459d0))
         @assertRelativelyEqual(3185.6731730807096d0,  soln%eq_soln(2)%T, tol)
         @assertRelativelyEqual(-6.1821344716750914d0, soln%eq_soln(2)%ln_nj(1), tol)
         @assertRelativelyEqual(-3.7758918668469459d0, soln%eq_soln(2)%ln_nj(2), tol)
@@ -72,6 +80,10 @@ contains
         @assertEqual(0.0d0, soln%eq_soln(1)%nj(11), tol)
 
         ! Test the exit conditions (pi_p)
+        ! write(*,*) "st3 T = ", soln%eq_soln(3)%T, " rel = ", &
+        !    abs((2567.3401804441942d0 - soln%eq_soln(3)%T)/2567.3401804441942d0)
+        ! write(*,*) "st3 ln_nj(H2) = ", soln%eq_soln(3)%ln_nj(2), " rel = ", &
+        !    abs((-3.7880413204534467d0 - soln%eq_soln(3)%ln_nj(2))/(-3.7880413204534467d0))
         @assertRelativelyEqual(2567.3401804441942d0,  soln%eq_soln(3)%T, tol)
         @assertRelativelyEqual(-7.4110702063645491d0, soln%eq_soln(3)%ln_nj(1), tol)
         @assertRelativelyEqual(-3.7880413204534467d0, soln%eq_soln(3)%ln_nj(2), tol)

source/shock.f90

@@ -145,6 +145,71 @@ subroutine ShockSolver_update_solution(self, soln, X1, X2, p21, t21, iter)
 
     end subroutine
 
+    logical function ShockSolver_state_valid(soln, idx) result(valid)
+        type(ShockSolution), intent(in) :: soln
+        integer, intent(in) :: idx
+
+        valid = .false.
+        if (idx < 1 .or. idx > soln%num_pts) return
+
+        valid = soln%eq_soln(idx)%converged .and. &
+                soln%eq_soln(idx)%T > 0.0d0 .and. &
+                soln%eq_soln(idx)%n > 0.0d0
+    end function
+
+    subroutine ShockSolver_fail_state(soln, idx, message)
+        type(ShockSolution), intent(inout) :: soln
+        integer, intent(in) :: idx
+        character(*), intent(in) :: message
+
+        if (len_trim(message) > 0) call log_warning(trim(message))
+        call log_warning("Shock calculation stopped after the last valid point.")
+
+        soln%converged = .false.
+        soln%pressure(idx) = 0.0d0
+        soln%mach(idx) = 0.0d0
+        soln%u(idx) = 0.0d0
+        soln%v_sonic(idx) = 0.0d0
+        soln%eq_soln(idx)%T = 0.0d0
+        soln%eq_soln(idx)%n = 0.0d0
+        soln%eq_soln(idx)%converged = .false.
+        soln%eq_partials(idx)%dlnV_dlnP = 0.0d0
+        soln%eq_partials(idx)%dlnV_dlnT = 0.0d0
+        soln%eq_partials(idx)%gamma_s = 0.0d0
+
+        select case (idx)
+            case (2)
+                soln%rho12 = 0.0d0
+                soln%p21 = 0.0d0
+                soln%t21 = 0.0d0
+                soln%M21 = 0.0d0
+                soln%v2 = 0.0d0
+                if (soln%num_pts >= 3) then
+                    soln%pressure(3) = 0.0d0
+                    soln%mach(3) = 0.0d0
+                    soln%u(3) = 0.0d0
+                    soln%v_sonic(3) = 0.0d0
+                    soln%eq_soln(3)%T = 0.0d0
+                    soln%eq_soln(3)%n = 0.0d0
+                    soln%eq_soln(3)%converged = .false.
+                    soln%eq_partials(3)%dlnV_dlnP = 0.0d0
+                    soln%eq_partials(3)%dlnV_dlnT = 0.0d0
+                    soln%eq_partials(3)%gamma_s = 0.0d0
+                    soln%rho52 = 0.0d0
+                    soln%p52 = 0.0d0
+                    soln%t52 = 0.0d0
+                    soln%M52 = 0.0d0
+                    soln%u5_p_v2 = 0.0d0
+                end if
+            case (3)
+                soln%rho52 = 0.0d0
+                soln%p52 = 0.0d0
+                soln%t52 = 0.0d0
+                soln%M52 = 0.0d0
+                soln%u5_p_v2 = 0.0d0
+        end select
+    end subroutine
+
     subroutine ShockSolver_solve_incident(self, soln, weights, T0, P0)
         ! Solve the incident shock conditions
 
@@ -198,6 +263,11 @@ subroutine ShockSolver_solve_incident(self, soln, weights, T0, P0)
 
         soln%pressure(idx) = p21*P0
         call self%eq_solver%solve(soln%eq_soln(idx), "hp", h0, soln%pressure(idx), weights, partials=soln%eq_partials(idx))
+        if (.not. ShockSolver_state_valid(soln, idx)) then
+            call ShockSolver_fail_state(soln, idx, &
+                                        "ShockSolver_solve_incident: incident equilibrium initialization failed.")
+            return
+        end if
 
         t21 = soln%eq_soln(idx)%T/T0
         ttmax = 1.05*T_gas_max/T0
@@ -209,6 +279,11 @@ subroutine ShockSolver_solve_incident(self, soln, weights, T0, P0)
             T2 = t21*T0
 
             call self%eq_solver%solve(soln%eq_soln(idx), "tp", T2, soln%pressure(idx), weights, partials=soln%eq_partials(idx))
+            if (.not. ShockSolver_state_valid(soln, idx)) then
+                call ShockSolver_fail_state(soln, idx, &
+                                            "ShockSolver_solve_incident: incident equilibrium iteration failed.")
+                return
+            end if
 
             ! Update properties after the equilibrium shock
             wm_k = 1.0d0/soln%eq_soln(idx)%n
@@ -244,10 +319,8 @@ subroutine ShockSolver_solve_incident(self, soln, weights, T0, P0)
             call self%update_solution(soln, X(1), X(2), p21, t21, i)
 
             if (i == 1 .and. .not. soln%converged .and. t21 >= ttmax) then
-                call log_warning("ShockSolver_solve_incident: first-iteration update hit " // &
-                                 "temperature cap; marking incident point as failed.")
-                soln%eq_soln(idx)%T = 0.0d0
-                soln%pressure(idx) = 0.0d0
+                call ShockSolver_fail_state(soln, idx, "ShockSolver_solve_incident: first-iteration update hit " // &
+                                            "temperature cap; marking incident point as failed.")
                 return
             end if
 
@@ -268,14 +341,7 @@ subroutine ShockSolver_solve_incident(self, soln, weights, T0, P0)
 
         ! Not converged; compute shock properties
         if (.not. soln%converged) then
-            soln%rho12 = rho12
-            soln%p21 = p21
-            soln%t21 = t21
-            soln%u(idx) = u1*rho12
-            soln%M21 = wm_k/wm
-            soln%mach(idx) = soln%M21*mach1
-            soln%v_sonic(idx) = (R*T2*(cp/(cp - 1.0/wm_k))/wm_k)**0.5d0
-            soln%v2 = u1 - soln%u(idx)
+            call ShockSolver_fail_state(soln, idx, "ShockSolver_solve_incident: incident equilibrium solve did not converge.")
         end if
 
     end subroutine
@@ -394,10 +460,8 @@ subroutine ShockSolver_solve_incident_frozen(self, soln, weights, T0, P0)
             call self%update_solution(soln, X(1), X(2), p21, t21, i)
 
             if (i == 1 .and. .not. soln%converged .and. t21 >= ttmax) then
-                call log_warning("ShockSolver_solve_incident_frozen: first-iteration update hit " // &
-                                 "temperature cap; marking incident point as failed.")
-                soln%eq_soln(idx)%T = 0.0d0
-                soln%pressure(idx) = 0.0d0
+                call ShockSolver_fail_state(soln, idx, "ShockSolver_solve_incident_frozen: first-iteration update hit " // &
+                                            "temperature cap; marking incident point as failed.")
                 return
             end if
 
@@ -419,14 +483,7 @@ subroutine ShockSolver_solve_incident_frozen(self, soln, weights, T0, P0)
 
         ! Not converged; compute shock properties
         if (.not. soln%converged) then
-            soln%rho12 = rho12
-            soln%p21 = p21
-            soln%t21 = t21
-            soln%u(idx) = u1*rho12
-            soln%M21 = wm_k/wm
-            soln%mach(idx) = soln%M21*mach1
-            soln%v_sonic(idx) = (R*T2*(cp/(cp - 1.0/wm_k))/wm_k)**0.5d0
-            soln%v2 = u1 - soln%u(idx)
+            call ShockSolver_fail_state(soln, idx, "ShockSolver_solve_incident_frozen: incident frozen solve did not converge.")
         end if
 
     end subroutine
@@ -456,7 +513,7 @@ subroutine ShockSolver_solve_reflected(self, soln, weights, T0, P0)
         real(dp) :: G(2, 3)               ! Solution matrix
         real(dp) :: X(3)                  ! Solution vector
         real(dp) :: dlnV_dlnP, dlnV_dlnT  ! Partial derivatives
-        real(dp) :: rho12                 ! Ratios of chemical potential and density across the incident shock
+        real(dp) :: rho12                 ! Incident-shock density ratio carried into reflected solve
         real(dp) :: mu25rt, rho52         ! Ratios of chemical potential and density across the reflected shock
         real(dp) :: tmp                   ! Intermediate variabls
         real(dp), allocatable :: nj_g(:)  ! Total/gas species concentrations [kmol-per-kg]
@@ -492,6 +549,11 @@ subroutine ShockSolver_solve_reflected(self, soln, weights, T0, P0)
             T5 = t52*T2
 
             call self%eq_solver%solve(soln%eq_soln(idx), "tp", T5, soln%pressure(idx), weights, partials=soln%eq_partials(idx))
+            if (.not. ShockSolver_state_valid(soln, idx)) then
+                call ShockSolver_fail_state(soln, idx, &
+                                            "ShockSolver_solve_reflected: reflected equilibrium iteration failed.")
+                return
+            end if
 
             ! Update properties after the equilibrium shock
             wm_k = 1.0d0/soln%eq_soln(idx)%n
@@ -531,10 +593,8 @@ subroutine ShockSolver_solve_reflected(self, soln, weights, T0, P0)
             call self%update_solution(soln, X(1), X(2), p52, t52, i)
 
             if (i == 1 .and. .not. soln%converged .and. t52 >= ttmax) then
-                call log_warning("ShockSolver_solve_reflected: first-iteration update hit " // &
-                                 "temperature cap; marking reflected point as failed.")
-                soln%eq_soln(idx)%T = 0.0d0
-                soln%pressure(idx) = 0.0d0
+                call ShockSolver_fail_state(soln, idx, "ShockSolver_solve_reflected: first-iteration update hit " // &
+                                            "temperature cap; marking reflected point as failed.")
                 return
             end if
 
@@ -555,14 +615,7 @@ subroutine ShockSolver_solve_reflected(self, soln, weights, T0, P0)
 
         ! Not converged; compute shock properties
         if (.not. soln%converged) then
-            soln%rho52 = rho52
-            soln%p52 = p52
-            soln%t52 = t52
-            soln%u(idx) = (u1 - u1*rho12)/rho52
-            soln%M52 = wm_k/wm
-            soln%mach(idx) = soln%M52*soln%mach(2)
-            soln%u5_p_v2 = (u1 - u1*rho12)*(1.0d0-1.0d0/rho52)
-            soln%v_sonic(idx) = (R*T5*(cp/(cp - 1.0/wm_k))/wm_k)**0.5d0
+            call ShockSolver_fail_state(soln, idx, "ShockSolver_solve_reflected: reflected equilibrium solve did not converge.")
         end if
 
     end subroutine
@@ -594,6 +647,7 @@ subroutine ShockSolver_solve_reflected_frozen(self, soln, weights, T0, P0)
         real(dp) :: dlnV_dlnP, dlnV_dlnT  ! Partial derivatives
         real(dp) :: rho12                 ! Ratios of chemical potential and density across the incident shock
         real(dp) :: mu25rt, rho52         ! Ratios of chemical potential and density across the reflected shock
+        real(dp) :: rho25_inv             ! Reflected-shock inverse density ratio used in the Newton system
         real(dp) :: tmp                   ! Intermediate variabls
         integer, parameter :: max_iter = 60
         real(dp), parameter :: T_gas_max = 20000.d0  ! Max gas temperature in the thermo database [K]
@@ -652,8 +706,8 @@ subroutine ShockSolver_solve_reflected_frozen(self, soln, weights, T0, P0)
             dlnV_dlnP = soln%eq_partials(idx)%dlnV_dlnP
             dlnV_dlnT = soln%eq_partials(idx)%dlnV_dlnT
 
-            rho12 = wm*t52/(wm_k*p52)
-            rho52 = 1./rho12
+            rho25_inv = wm*t52/(wm_k*p52)
+            rho52 = 1./rho25_inv
             soln%rho52 = rho52
             tmp = -mu25rt*rho52/(rho52 - 1.0d0)**2
 
@@ -677,10 +731,8 @@ subroutine ShockSolver_solve_reflected_frozen(self, soln, weights, T0, P0)
             call self%update_solution(soln, X(1), X(2), p52, t52, i)
 
             if (i == 1 .and. .not. soln%converged .and. t52 >= ttmax) then
-                call log_warning("ShockSolver_solve_reflected_frozen: first-iteration update hit " // &
-                                 "temperature cap; marking reflected point as failed.")
-                soln%eq_soln(idx)%T = 0.0d0
-                soln%pressure(idx) = 0.0d0
+                call ShockSolver_fail_state(soln, idx, "ShockSolver_solve_reflected_frozen: first-iteration update hit " // &
+                                            "temperature cap; marking reflected point as failed.")
                 return
             end if
 
@@ -701,14 +753,7 @@ subroutine ShockSolver_solve_reflected_frozen(self, soln, weights, T0, P0)
 
         ! Not converged; compute shock properties
         if (.not. soln%converged) then
-            soln%rho52 = rho52
-            soln%p52 = p52
-            soln%t52 = t52
-            soln%u(idx) = soln%v2/rho52
-            soln%M52 = wm_k/wm
-            soln%mach(idx) = soln%M52*soln%mach(2)
-            soln%u5_p_v2 = (u1 - u1*rho12)*(1.0d0-1.0d0/rho52)
-            soln%v_sonic(idx) = (R*T5*(cp/(cp - 1.0/wm_k))/wm_k)**0.5d0
+            call ShockSolver_fail_state(soln, idx, "ShockSolver_solve_reflected_frozen: reflected frozen solve did not converge.")
         end if
 
     end subroutine