From 967c1284221c54d863c69640350b2b000545c5e8 Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Sun, 29 Jun 2025 02:56:17 -0400 Subject: [PATCH 1/9] s_compute_wave_speed --- src/common/m_variables_conversion.fpp | 115 ++++++++++++++ src/simulation/m_riemann_solvers.fpp | 218 +++----------------------- 2 files changed, 135 insertions(+), 198 deletions(-) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 949eac92c..c150feaf2 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -44,6 +44,7 @@ module m_variables_conversion #ifndef MFC_PRE_PROCESS s_compute_speed_of_sound, & s_compute_fast_magnetosonic_speed, & + s_compute_wave_speed, & #endif s_finalize_variables_conversion_module @@ -1709,4 +1710,118 @@ contains end subroutine s_compute_fast_magnetosonic_speed #endif +#ifndef MFC_PRE_PROCESS + subroutine s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, & + c_L, c_R, c_avg, c_fast_L, c_fast_R, G_L, G_R, & + tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & + s_L, s_R, s_S, s_M, s_P, idx, idx_tau) + + ! Computes the wave speeds for the Riemann solver +#ifdef _CRAYFTN + !DIR$ INLINEALWAYS s_compute_wave_speed +#else + !$acc routine seq +#endif + + ! Input parameters + integer, intent(in) :: wave_speeds + integer, intent(in) :: idx, idx_tau + real(wp), intent(in) :: rho_L, rho_R + real(wp), dimension(:), intent(in) :: vel_L, vel_R, tau_e_L, tau_e_R + real(wp), intent(in) :: pres_L, pres_R, c_L, c_R + real(wp), intent(in) :: gamma_L, gamma_R, pi_inf_L, pi_inf_R + real(wp), intent(in) :: rho_avg, c_avg + real(wp), intent(in) :: c_fast_L, c_fast_R + real(wp), intent(in) :: G_L, G_R + + ! Local variables + real(wp) :: pres_SL, pres_SR, Ms_L, Ms_R + + ! Output parameters + real(wp), intent(out) :: s_L, s_R, s_S, s_M, s_P + + if (wave_speeds == 1) then + if (elasticity) then + s_L = min(vel_L(idx) - sqrt(c_L*c_L + & + (((4_wp*G_L)/3_wp) + tau_e_L(idx_tau))/rho_L), vel_R(idx) - sqrt(c_R*c_R + & + (((4_wp*G_R)/3_wp) + tau_e_R(idx_tau))/rho_R)) + s_R = max(vel_R(idx) + sqrt(c_R*c_R + & + (((4_wp*G_R)/3_wp) + tau_e_R(idx_tau))/rho_R), vel_L(idx) + sqrt(c_L*c_L + & + (((4_wp*G_L)/3_wp) + tau_e_L(idx_tau))/rho_L)) + s_S = (pres_R - tau_e_R(idx_tau) - pres_L + & + tau_e_L(idx_tau) + rho_L*vel_L(idx)*(s_L - vel_L(idx)) - & + rho_R*vel_R(idx)*(s_R - vel_R(idx)))/(rho_L*(s_L - vel_L(idx)) - & + rho_R*(s_R - vel_R(idx))) + else if (mhd) then + s_L = min(vel_L(idx) - c_fast_L, vel_R(idx) - c_fast_R) + s_R = max(vel_R(idx) + c_fast_R, vel_L(idx) + c_fast_L) + s_S = (pres_R - pres_L + rho_L*vel_L(idx)* & + (s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) & + /(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx))) + else if (hypoelasticity) then + s_L = min(vel_L(idx) - sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + & + tau_e_L(idx_tau))/rho_L) & + , vel_R(idx) - sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) + & + tau_e_R(idx_tau))/rho_R)) + s_R = max(vel_R(idx) + sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) + & + tau_e_R(idx_tau))/rho_R) & + , vel_L(idx) + sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + & + tau_e_L(idx_tau))/rho_L)) + s_S = (pres_R - pres_L + rho_L*vel_L(idx)* & + (s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) & + /(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx))) + else if (hyperelasticity) then + s_L = min(vel_L(idx) - sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L) & + , vel_R(idx) - sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R)) + s_R = max(vel_R(idx) + sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R) & + , vel_L(idx) + sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L)) + s_S = (pres_R - pres_L + rho_L*vel_L(idx)* & + (s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) & + /(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx))) + else + s_L = min(vel_L(idx) - c_L, vel_R(idx) - c_R) + s_R = max(vel_R(idx) + c_R, vel_L(idx) + c_L) + s_S = (pres_R - pres_L + rho_L*vel_L(idx)* & + (s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) & + /(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx))) + end if + else if (wave_speeds == 2) then + pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg*(vel_L(idx) - vel_R(idx))) + pres_SR = pres_SL + Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + (pres_SL/pres_L - 1._wp)*pres_L/ & + ((pres_L + pi_inf_L/(1._wp + gamma_L))))) + Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + (pres_SR/pres_R - 1._wp)*pres_R/ & + ((pres_R + pi_inf_R/(1._wp + gamma_R))))) + s_L = vel_L(idx) - c_L*Ms_L + s_R = vel_R(idx) + c_R*Ms_R + s_S = 5e-1_wp*((vel_L(idx) + vel_R(idx)) + (pres_L - pres_R)/(rho_avg*c_avg)) + end if + + ! ! follows Einfeldt et al. + ! s_M/P = min/max(0.,s_L/R) + s_M = min(0._wp, s_L) + s_P = max(0._wp, s_R) + +#ifdef DEBUG + ! Check for potential issues in wave speed calculation + if (s_R <= s_L) then + print *, 'WARNING: Wave speed issue detected in s_compute_wave_speed' + print *, 'Left wave speed >= Right wave speed:', s_L, s_R + print *, 'Input velocities :', vel_L(idx), vel_R(idx) + print *, 'Sound speeds:', c_L, c_R + print *, 'Densities:', rho_L, rho_R + print *, 'Pressures:', pres_L, pres_R + print *, 'Wave speeds method:', wave_speeds + if (elasticity .or. hypoelasticity .or. hyperelasticity) then + print *, 'Shear moduli:', G_L, G_R + end if + call s_mpi_abort('Error: Invalid wave speeds in s_compute_wave_speed') + end if +#endif + + end subroutine s_compute_wave_speed +#endif + end module m_variables_conversion diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 9c65e51e2..dbe96f7ff 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -682,60 +682,10 @@ contains end do end if - if (wave_speeds == 1) then - if (mhd) then - s_L = min(vel_L(dir_idx(1)) - c_fast%L, vel_R(dir_idx(1)) - c_fast%R) - s_R = max(vel_R(dir_idx(1)) + c_fast%R, vel_L(dir_idx(1)) + c_fast%L) - elseif (hypoelasticity) then - s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + & - (((4._wp*G_L)/3._wp) + & - tau_e_L(dir_idx_tau(1)))/rho_L) & - , vel_R(dir_idx(1)) - sqrt(c_R*c_R + & - (((4._wp*G_R)/3._wp) + & - tau_e_R(dir_idx_tau(1)))/rho_R)) - s_R = max(vel_R(dir_idx(1)) + sqrt(c_R*c_R + & - (((4._wp*G_R)/3._wp) + & - tau_e_R(dir_idx_tau(1)))/rho_R) & - , vel_L(dir_idx(1)) + sqrt(c_L*c_L + & - (((4._wp*G_L)/3._wp) + & - tau_e_L(dir_idx_tau(1)))/rho_L)) - else if (hyperelasticity) then - s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L) & - , vel_R(dir_idx(1)) - sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R)) - s_R = max(vel_R(dir_idx(1)) + sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R) & - , vel_L(dir_idx(1)) + sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L)) - else - s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) - s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) - end if - - s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))* & - (s_L - vel_L(dir_idx(1))) - & - rho_R*vel_R(dir_idx(1))* & - (s_R - vel_R(dir_idx(1)))) & - /(rho_L*(s_L - vel_L(dir_idx(1))) - & - rho_R*(s_R - vel_R(dir_idx(1)))) - elseif (wave_speeds == 2) then - pres_SL = 5.e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) - - pres_SR = pres_SL - - Ms_L = max(1._wp, sqrt(1._wp + ((5.e-1_wp + gamma_L)/(1._wp + gamma_L))* & - (pres_SL/pres_L - 1._wp)*pres_L/ & - ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5.e-1_wp + gamma_R)/(1._wp + gamma_R))* & - (pres_SR/pres_R - 1._wp)*pres_R/ & - ((pres_R + pi_inf_R/(1._wp + gamma_R))))) - - s_L = vel_L(dir_idx(1)) - c_L*Ms_L - s_R = vel_R(dir_idx(1)) + c_R*Ms_R - - s_S = 5.e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) - end if + call s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, & + c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, & + tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & + s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) @@ -1488,47 +1438,10 @@ contains end if ! COMPUTING THE DIRECT WAVE SPEEDS - if (wave_speeds == 1) then - if (elasticity) then - s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + & - (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L), vel_R(dir_idx(1)) - sqrt(c_R*c_R + & - (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)) - s_R = max(vel_R(dir_idx(1)) + sqrt(c_R*c_R + & - (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R), vel_L(dir_idx(1)) + sqrt(c_L*c_L + & - (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)) - s_S = (pres_R - tau_e_R(dir_idx_tau(1)) - pres_L + & - tau_e_L(dir_idx_tau(1)) + rho_L*vel_L(idx1)*(s_L - vel_L(idx1)) - & - rho_R*vel_R(idx1)*(s_R - vel_R(idx1)))/(rho_L*(s_L - vel_L(idx1)) - & - rho_R*(s_R - vel_R(idx1))) - else - s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) - s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) - s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))* & - (s_L - vel_L(dir_idx(1))) - rho_R*vel_R(dir_idx(1))*(s_R - vel_R(dir_idx(1)))) & - /(rho_L*(s_L - vel_L(dir_idx(1))) - rho_R*(s_R - vel_R(dir_idx(1)))) - - end if - elseif (wave_speeds == 2) then - pres_SL = 5.e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) - - pres_SR = pres_SL - - Ms_L = max(1._wp, sqrt(1._wp + ((5.e-1_wp + gamma_L)/(1._wp + gamma_L))* & - (pres_SL/pres_L - 1._wp)*pres_L/ & - ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5.e-1_wp + gamma_R)/(1._wp + gamma_R))* & - (pres_SR/pres_R - 1._wp)*pres_R/ & - ((pres_R + pi_inf_R/(1._wp + gamma_R))))) - - s_L = vel_L(dir_idx(1)) - c_L*Ms_L - s_R = vel_R(dir_idx(1)) + c_R*Ms_R - - s_S = 5.e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) - end if + call s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, & + c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, & + tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & + s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) ! follows Einfeldt et al. ! s_M/P = min/max(0.,s_L/R) @@ -1798,37 +1711,10 @@ contains call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & vel_avg_rms, 0._wp, c_avg) - if (wave_speeds == 1) then - s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) - s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) - - s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))* & - (s_L - vel_L(dir_idx(1))) - & - rho_R*vel_R(dir_idx(1))* & - (s_R - vel_R(dir_idx(1)))) & - /(rho_L*(s_L - vel_L(dir_idx(1))) - & - rho_R*(s_R - vel_R(dir_idx(1)))) - elseif (wave_speeds == 2) then - pres_SL = 5.e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) - - pres_SR = pres_SL - - Ms_L = max(1._wp, sqrt(1._wp + ((5.e-1_wp + gamma_L)/(1._wp + gamma_L))* & - (pres_SL/pres_L - 1._wp)*pres_L/ & - ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5.e-1_wp + gamma_R)/(1._wp + gamma_R))* & - (pres_SR/pres_R - 1._wp)*pres_R/ & - ((pres_R + pi_inf_R/(1._wp + gamma_R))))) - - s_L = vel_L(dir_idx(1)) - c_L*Ms_L - s_R = vel_R(dir_idx(1)) + c_R*Ms_R - - s_S = 5.e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) - end if + call s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, & + c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, & + tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & + s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) ! follows Einfeldt et al. ! s_M/P = min/max(0.,s_L/R) @@ -2227,37 +2113,10 @@ contains @:compute_low_Mach_correction() end if - if (wave_speeds == 1) then - s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) - s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) - - s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))* & - (s_L - vel_L(dir_idx(1))) - & - rho_R*vel_R(dir_idx(1))* & - (s_R - vel_R(dir_idx(1)))) & - /(rho_L*(s_L - vel_L(dir_idx(1))) - & - rho_R*(s_R - vel_R(dir_idx(1)))) - elseif (wave_speeds == 2) then - pres_SL = 5.e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) - - pres_SR = pres_SL - - Ms_L = max(1._wp, sqrt(1._wp + ((5.e-1_wp + gamma_L)/(1._wp + gamma_L))* & - (pres_SL/pres_L - 1._wp)*pres_L/ & - ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5.e-1_wp + gamma_R)/(1._wp + gamma_R))* & - (pres_SR/pres_R - 1._wp)*pres_R/ & - ((pres_R + pi_inf_R/(1._wp + gamma_R))))) - - s_L = vel_L(dir_idx(1)) - c_L*Ms_L - s_R = vel_R(dir_idx(1)) + c_R*Ms_R - - s_S = 5.e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) - end if + call s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, & + c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, & + tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & + s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) ! follows Einfeldt et al. ! s_M/P = min/max(0.,s_L/R) @@ -2696,47 +2555,10 @@ contains @:compute_low_Mach_correction() end if - if (wave_speeds == 1) then - if (elasticity) then - s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + & - (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L), vel_R(dir_idx(1)) - sqrt(c_R*c_R + & - (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)) - s_R = max(vel_R(dir_idx(1)) + sqrt(c_R*c_R + & - (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R), vel_L(dir_idx(1)) + sqrt(c_L*c_L + & - (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)) - s_S = (pres_R - tau_e_R(dir_idx_tau(1)) - pres_L + & - tau_e_L(dir_idx_tau(1)) + rho_L*vel_L(idx1)*(s_L - vel_L(idx1)) - & - rho_R*vel_R(idx1)*(s_R - vel_R(idx1)))/(rho_L*(s_L - vel_L(idx1)) - & - rho_R*(s_R - vel_R(idx1))) - else - s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) - s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) - s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))* & - (s_L - vel_L(dir_idx(1))) - rho_R*vel_R(dir_idx(1))*(s_R - vel_R(dir_idx(1)))) & - /(rho_L*(s_L - vel_L(dir_idx(1))) - rho_R*(s_R - vel_R(dir_idx(1)))) - - end if - elseif (wave_speeds == 2) then - pres_SL = 5.e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(idx1) - & - vel_R(idx1))) - - pres_SR = pres_SL - - Ms_L = max(1._wp, sqrt(1._wp + ((5.e-1_wp + gamma_L)/(1._wp + gamma_L))* & - (pres_SL/pres_L - 1._wp)*pres_L/ & - ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5.e-1_wp + gamma_R)/(1._wp + gamma_R))* & - (pres_SR/pres_R - 1._wp)*pres_R/ & - ((pres_R + pi_inf_R/(1._wp + gamma_R))))) - - s_L = vel_L(idx1) - c_L*Ms_L - s_R = vel_R(idx1) + c_R*Ms_R - - s_S = 5.e-1_wp*((vel_L(idx1) + vel_R(idx1)) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) - end if + call s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, & + c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, & + tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & + s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) ! follows Einfeldt et al. ! s_M/P = min/max(0.,s_L/R) From 89280ad2dc814e7d488740c8241f8316418a446f Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Sun, 29 Jun 2025 03:12:59 -0400 Subject: [PATCH 2/9] low-stakes refactoring & lint/format --- src/simulation/m_riemann_solvers.fpp | 632 +++++++-------------------- 1 file changed, 154 insertions(+), 478 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index dbe96f7ff..8442dad6e 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -683,9 +683,9 @@ contains end if call s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, & - c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, & - tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & - s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) + c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, & + tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & + s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) @@ -3269,13 +3269,20 @@ contains dqR_prim_dz_vf, & norm_dir, ix, iy, iz) - real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf + real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), target, intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf + real(wp), dimension(:, :, :, :), pointer :: qL_prim_rs_vf, qR_prim_rs_vf type(scalar_field), & allocatable, dimension(:), & - intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, & - dqL_prim_dy_vf, dqR_prim_dy_vf, & - dqL_prim_dz_vf, dqR_prim_dz_vf + target, intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, & + dqL_prim_dy_vf, dqR_prim_dy_vf, & + dqL_prim_dz_vf, dqR_prim_dz_vf, & + qL_prim_vf, qR_prim_vf + type(scalar_field), & + dimension(:), & + pointer :: dqL_prim_d_vf, dqR_prim_d_vf + + integer :: end_val, bc_beg, bc_end integer, intent(in) :: norm_dir type(int_bounds_info), intent(in) :: ix, iy, iz @@ -3285,12 +3292,30 @@ contains if (norm_dir == 1) then is1 = ix; is2 = iy; is3 = iz dir_idx = (/1, 2, 3/); dir_flg = (/1._wp, 0._wp, 0._wp/) - elseif (norm_dir == 2) then + bc_beg = bc_x%beg; bc_end = bc_x%end + end_val = m + qL_prim_rs_vf => qL_prim_rsx_vf + qR_prim_rs_vf => qR_prim_rsx_vf + dqL_prim_d_vf => dqL_prim_dx_vf + dqR_prim_d_vf => dqR_prim_dx_vf + else if (norm_dir == 2) then is1 = iy; is2 = ix; is3 = iz dir_idx = (/2, 1, 3/); dir_flg = (/0._wp, 1._wp, 0._wp/) + bc_beg = bc_y%beg; bc_end = bc_y%end + end_val = n + qL_prim_rs_vf => qL_prim_rsy_vf + qR_prim_rs_vf => qR_prim_rsy_vf + dqL_prim_d_vf => dqL_prim_dy_vf + dqR_prim_d_vf => dqR_prim_dy_vf else is1 = iz; is2 = iy; is3 = ix dir_idx = (/3, 1, 2/); dir_flg = (/0._wp, 0._wp, 1._wp/) + bc_beg = bc_z%beg; bc_end = bc_z%end + end_val = p + qL_prim_rs_vf => qL_prim_rsz_vf + qR_prim_rs_vf => qR_prim_rsz_vf + dqL_prim_d_vf => dqL_prim_dz_vf + dqR_prim_d_vf => dqR_prim_dz_vf end if !$acc update device(is1, is2, is3) @@ -3309,315 +3334,83 @@ contains !$acc update device(isx, isy, isz) ! for stuff in the same module !$acc update device(dir_idx, dir_flg, dir_idx_tau) ! for stuff in different modules - ! Population of Buffers in x-direction - if (norm_dir == 1) then - - if (bc_x%beg == BC_RIEMANN_EXTRAP) then ! Riemann state extrap. BC at beginning - !$acc parallel loop collapse(3) gang vector default(present) - do i = 1, sys_size - do l = is3%beg, is3%end - do k = is2%beg, is2%end - qL_prim_rsx_vf(-1, k, l, i) = & - qR_prim_rsx_vf(0, k, l, i) - end do - end do - end do - - if (viscous) then - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do k = isy%beg, isy%end - - dqL_prim_dx_vf(i)%sf(-1, k, l) = & - dqR_prim_dx_vf(i)%sf(0, k, l) - end do - end do - end do - - if (n > 0) then - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do k = isy%beg, isy%end - - dqL_prim_dy_vf(i)%sf(-1, k, l) = & - dqR_prim_dy_vf(i)%sf(0, k, l) - end do - end do - end do - - if (p > 0) then - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do k = isy%beg, isy%end - - dqL_prim_dz_vf(i)%sf(-1, k, l) = & - dqR_prim_dz_vf(i)%sf(0, k, l) - end do - end do - end do - end if - - end if - - end if - - end if - - if (bc_x%end == BC_RIEMANN_EXTRAP) then ! Riemann state extrap. BC at end - - !$acc parallel loop collapse(3) gang vector default(present) - do i = 1, sys_size - do l = is3%beg, is3%end - do k = is2%beg, is2%end - qR_prim_rsx_vf(m + 1, k, l, i) = & - qL_prim_rsx_vf(m, k, l, i) - end do + ! Population of Buffers in x/y/z-direction + if (bc_beg == BC_RIEMANN_EXTRAP) then ! Riemann state extrap. BC at beginning + !$acc parallel loop collapse(3) gang vector default(present) + do i = 1, sys_size + do l = is3%beg, is3%end + do k = is2%beg, is2%end + qL_prim_rs_vf(-1, k, l, i) = qR_prim_rs_vf(0, k, l, i) end do end do - - if (viscous) then - - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do k = isy%beg, isy%end - - dqR_prim_dx_vf(i)%sf(m + 1, k, l) = & - dqL_prim_dx_vf(i)%sf(m, k, l) - end do - end do - end do - - if (n > 0) then - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do k = isy%beg, isy%end - - dqR_prim_dy_vf(i)%sf(m + 1, k, l) = & - dqL_prim_dy_vf(i)%sf(m, k, l) - end do - end do - end do - - if (p > 0) then - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do k = isy%beg, isy%end - - dqR_prim_dz_vf(i)%sf(m + 1, k, l) = & - dqL_prim_dz_vf(i)%sf(m, k, l) - end do - end do - end do - end if - - end if - - end if - - end if - ! END: Population of Buffers in x-direction - - ! Population of Buffers in y-direction - elseif (norm_dir == 2) then - - if (bc_y%beg == BC_RIEMANN_EXTRAP) then ! Riemann state extrap. BC at beginning + end do + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) - do i = 1, sys_size - do l = is3%beg, is3%end - do k = is2%beg, is2%end - qL_prim_rsy_vf(-1, k, l, i) = & - qR_prim_rsy_vf(0, k, l, i) + do i = momxb, momxe + do l = isz%beg, isz%end + do k = isy%beg, isy%end + if (norm_dir == 1) then + dqL_prim_dx_vf(i)%sf(-1, k, l) = dqR_prim_dx_vf(i)%sf(0, k, l) + if (n > 0) then + dqL_prim_dy_vf(i)%sf(-1, k, l) = dqR_prim_dy_vf(i)%sf(0, k, l) + if (p > 0) then + dqL_prim_dz_vf(i)%sf(-1, k, l) = dqR_prim_dz_vf(i)%sf(0, k, l) + end if + end if + else if (norm_dir == 2) then + dqL_prim_dx_vf(i)%sf(j, -1, l) = dqR_prim_dx_vf(i)%sf(j, 0, l) + dqL_prim_dy_vf(i)%sf(j, -1, l) = dqR_prim_dy_vf(i)%sf(j, 0, l) + if (p > 0) then + dqL_prim_dz_vf(i)%sf(j, -1, l) = dqR_prim_dz_vf(i)%sf(j, 0, l) + end if + else + dqL_prim_dx_vf(i)%sf(j, k, -1) = dqR_prim_dx_vf(i)%sf(j, k, 0) + dqL_prim_dy_vf(i)%sf(j, k, -1) = dqR_prim_dy_vf(i)%sf(j, k, 0) + dqL_prim_dz_vf(i)%sf(j, k, -1) = dqR_prim_dz_vf(i)%sf(j, k, 0) + end if end do end do end do - - if (viscous) then - - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do j = isx%beg, isx%end - dqL_prim_dx_vf(i)%sf(j, -1, l) = & - dqR_prim_dx_vf(i)%sf(j, 0, l) - end do - end do - end do - - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do j = isx%beg, isx%end - dqL_prim_dy_vf(i)%sf(j, -1, l) = & - dqR_prim_dy_vf(i)%sf(j, 0, l) - end do - end do - end do - - if (p > 0) then - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do j = isx%beg, isx%end - dqL_prim_dz_vf(i)%sf(j, -1, l) = & - dqR_prim_dz_vf(i)%sf(j, 0, l) - end do - end do - end do - end if - - end if - end if + end if - if (bc_y%end == BC_RIEMANN_EXTRAP) then ! Riemann state extrap. BC at end - - !$acc parallel loop collapse(3) gang vector default(present) - do i = 1, sys_size - do l = is3%beg, is3%end - do k = is2%beg, is2%end - qR_prim_rsy_vf(n + 1, k, l, i) = & - qL_prim_rsy_vf(n, k, l, i) - end do + if (bc_end == BC_RIEMANN_EXTRAP) then ! Riemann state extrap. BC at end + !$acc parallel loop collapse(3) gang vector default(present) + do i = 1, sys_size + do l = is3%beg, is3%end + do k = is2%beg, is2%end + qR_prim_rs_vf(end_val + 1, k, l, i) = qL_prim_rs_vf(end_val, k, l, i) end do end do - - if (viscous) then - - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do j = isx%beg, isx%end - dqR_prim_dx_vf(i)%sf(j, n + 1, l) = & - dqL_prim_dx_vf(i)%sf(j, n, l) - end do - end do - end do - - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do j = isx%beg, isx%end - dqR_prim_dy_vf(i)%sf(j, n + 1, l) = & - dqL_prim_dy_vf(i)%sf(j, n, l) - end do - end do - end do - - if (p > 0) then - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do l = isz%beg, isz%end - do j = isx%beg, isx%end - dqR_prim_dz_vf(i)%sf(j, n + 1, l) = & - dqL_prim_dz_vf(i)%sf(j, n, l) - end do - end do - end do - end if - - end if - - end if - ! END: Population of Buffers in y-direction - - ! Population of Buffers in z-direction - else - - if (bc_z%beg == BC_RIEMANN_EXTRAP) then ! Riemann state extrap. BC at beginning + end do + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) - do i = 1, sys_size - do l = is3%beg, is3%end - do k = is2%beg, is2%end - qL_prim_rsz_vf(-1, k, l, i) = & - qR_prim_rsz_vf(0, k, l, i) - end do - end do - end do - - if (viscous) then - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do k = isy%beg, isy%end - do j = isx%beg, isx%end - dqL_prim_dx_vf(i)%sf(j, k, -1) = & - dqR_prim_dx_vf(i)%sf(j, k, 0) - end do - end do - end do - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do k = isy%beg, isy%end - do j = isx%beg, isx%end - dqL_prim_dy_vf(i)%sf(j, k, -1) = & - dqR_prim_dy_vf(i)%sf(j, k, 0) - end do - end do - end do - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe + do i = momxb, momxe + do l = isz%beg, isz%end do k = isy%beg, isy%end - do j = isx%beg, isx%end - dqL_prim_dz_vf(i)%sf(j, k, -1) = & - dqR_prim_dz_vf(i)%sf(j, k, 0) - end do - end do - end do - end if - - end if - - if (bc_z%end == BC_RIEMANN_EXTRAP) then ! Riemann state extrap. BC at end - - !$acc parallel loop collapse(3) gang vector default(present) - do i = 1, sys_size - do l = is3%beg, is3%end - do k = is2%beg, is2%end - qR_prim_rsz_vf(p + 1, k, l, i) = & - qL_prim_rsz_vf(p, k, l, i) + if (norm_dir == 1) then + dqR_prim_dx_vf(i)%sf(end_val + 1, k, l) = dqL_prim_dx_vf(i)%sf(end_val, k, l) + if (n > 0) then + dqR_prim_dy_vf(i)%sf(end_val + 1, k, l) = dqL_prim_dy_vf(i)%sf(end_val, k, l) + if (p > 0) then + dqR_prim_dz_vf(i)%sf(end_val + 1, k, l) = dqL_prim_dz_vf(i)%sf(end_val, k, l) + end if + end if + else if (norm_dir == 2) then + dqR_prim_dx_vf(i)%sf(j, end_val + 1, l) = dqL_prim_dx_vf(i)%sf(j, end_val, l) + dqR_prim_dy_vf(i)%sf(j, end_val + 1, l) = dqL_prim_dy_vf(i)%sf(j, end_val, l) + if (p > 0) then + dqR_prim_dz_vf(i)%sf(j, end_val + 1, l) = dqL_prim_dz_vf(i)%sf(j, end_val, l) + end if + else + dqR_prim_dx_vf(i)%sf(j, k, end_val + 1) = dqL_prim_dx_vf(i)%sf(j, k, end_val) + dqR_prim_dy_vf(i)%sf(j, k, end_val + 1) = dqL_prim_dy_vf(i)%sf(j, k, end_val) + dqR_prim_dz_vf(i)%sf(j, k, end_val + 1) = dqL_prim_dz_vf(i)%sf(j, k, end_val) + end if end do end do end do - - if (viscous) then - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do k = isy%beg, isy%end - do j = isx%beg, isx%end - dqR_prim_dx_vf(i)%sf(j, k, p + 1) = & - dqL_prim_dx_vf(i)%sf(j, k, p) - end do - end do - end do - - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do k = isy%beg, isy%end - do j = isx%beg, isx%end - dqR_prim_dy_vf(i)%sf(j, k, p + 1) = & - dqL_prim_dy_vf(i)%sf(j, k, p) - end do - end do - end do - - !$acc parallel loop collapse(3) gang vector default(present) - do i = momxb, momxe - do k = isy%beg, isy%end - do j = isx%beg, isx%end - dqR_prim_dz_vf(i)%sf(j, k, p + 1) = & - dqL_prim_dz_vf(i)%sf(j, k, p) - end do - end do - end do - end if - end if - end if ! END: Population of Buffers in z-direction @@ -3653,94 +3446,42 @@ contains ! Reshaping Inputted Data in x-direction - if (norm_dir == 1) then - - if (viscous .or. (surface_tension)) then - - !$acc parallel loop collapse(4) gang vector default(present) - do i = momxb, E_idx - do l = is3%beg, is3%end - do k = is2%beg, is2%end - do j = is1%beg, is1%end - flux_src_vf(i)%sf(j, k, l) = 0._wp - end do - end do - end do - end do - end if - - if (qbmm) then - - !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, 4 - do l = is3%beg, is3%end - do k = is2%beg, is2%end - do j = is1%beg, is1%end + 1 - mom_sp_rsx_vf(j, k, l, i) = mom_sp(i)%sf(j, k, l) - end do - end do - end do - end do - end if - - ! Reshaping Inputted Data in y-direction - elseif (norm_dir == 2) then - - if (viscous .or. (surface_tension)) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = momxb, E_idx - do l = is3%beg, is3%end + if (viscous .or. (surface_tension)) then + !$acc parallel loop collapse(4) gang vector default(present) + do i = momxb, E_idx + do l = is3%beg, is3%end + do k = is2%beg, is2%end do j = is1%beg, is1%end - do k = is2%beg, is2%end + if (norm_dir == 1) then + flux_src_vf(i)%sf(j, k, l) = 0._wp + else if (norm_dir == 2) then flux_src_vf(i)%sf(k, j, l) = 0._wp - end do - end do - end do - end do - end if - - if (qbmm) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, 4 - do l = is3%beg, is3%end - do k = is2%beg, is2%end - do j = is1%beg, is1%end + 1 - mom_sp_rsy_vf(j, k, l, i) = mom_sp(i)%sf(k, j, l) - end do - end do - end do - end do - end if - - ! Reshaping Inputted Data in z-direction - else - - if (viscous .or. (surface_tension)) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = momxb, E_idx - do j = is1%beg, is1%end - do k = is2%beg, is2%end - do l = is3%beg, is3%end + else if (norm_dir == 3) then flux_src_vf(i)%sf(l, k, j) = 0._wp - end do + end if end do end do end do - end if + end do + end if - if (qbmm) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, 4 - do l = is3%beg, is3%end - do k = is2%beg, is2%end - do j = is1%beg, is1%end + 1 + if (qbmm) then + !$acc parallel loop collapse(4) gang vector default(present) + do i = 1, 4 + do l = is3%beg, is3%end + do k = is2%beg, is2%end + do j = is1%beg, is1%end + 1 + if (norm_dir == 1) then + mom_sp_rsx_vf(j, k, l, i) = mom_sp(i)%sf(j, k, l) + else if (norm_dir == 2) then + mom_sp_rsy_vf(j, k, l, i) = mom_sp(i)%sf(k, j, l) + else if (norm_dir == 3) then mom_sp_rsz_vf(j, k, l, i) = mom_sp(i)%sf(l, k, j) - end do + end if end do end do end do - end if - + end do end if end subroutine s_initialize_riemann_solver @@ -4127,145 +3868,80 @@ contains ! Reshaping Outputted Data in y-direction if (norm_dir == 2) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, sys_size - do l = is3%beg, is3%end - do j = is1%beg, is1%end - do k = is2%beg, is2%end - flux_vf(i)%sf(k, j, l) = & - flux_rsy_vf(j, k, l, i) - end do - end do - end do - end do - - if (cyl_coord) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, sys_size - do l = is3%beg, is3%end - do j = is1%beg, is1%end - do k = is2%beg, is2%end - flux_gsrc_vf(i)%sf(k, j, l) = & - flux_gsrc_rsy_vf(j, k, l, i) - end do - end do - end do - end do - end if - !$acc parallel loop collapse(3) gang vector default(present) do l = is3%beg, is3%end do j = is1%beg, is1%end do k = is2%beg, is2%end flux_src_vf(advxb)%sf(k, j, l) = & flux_src_rsy_vf(j, k, l, advxb) - end do - end do - end do - - if (riemann_solver == 1 .or. riemann_solver == 4) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = advxb + 1, advxe - do l = is3%beg, is3%end - do j = is1%beg, is1%end - do k = is2%beg, is2%end - flux_src_vf(i)%sf(k, j, l) = & - flux_src_rsy_vf(j, k, l, i) - end do + do i = 1, sys_size + flux_vf(i)%sf(k, j, l) = & + flux_rsy_vf(j, k, l, i) + if (cyl_coord) then + flux_gsrc_vf(i)%sf(k, j, l) = & + flux_gsrc_rsy_vf(j, k, l, i) + end if end do end do end do + end do - end if ! Reshaping Outputted Data in z-direction elseif (norm_dir == 3) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, sys_size - do j = is1%beg, is1%end - do k = is2%beg, is2%end - do l = is3%beg, is3%end - - flux_vf(i)%sf(l, k, j) = & - flux_rsz_vf(j, k, l, i) - end do - end do - end do - end do - if (grid_geometry == 3) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, sys_size - do j = is1%beg, is1%end - do k = is2%beg, is2%end - do l = is3%beg, is3%end - - flux_gsrc_vf(i)%sf(l, k, j) = & - flux_gsrc_rsz_vf(j, k, l, i) - end do - end do - end do - end do - end if - !$acc parallel loop collapse(3) gang vector default(present) do j = is1%beg, is1%end do k = is2%beg, is2%end do l = is3%beg, is3%end flux_src_vf(advxb)%sf(l, k, j) = & flux_src_rsz_vf(j, k, l, advxb) - end do - end do - end do - - if (riemann_solver == 1 .or. riemann_solver == 4) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = advxb + 1, advxe - do j = is1%beg, is1%end - do k = is2%beg, is2%end - do l = is3%beg, is3%end - flux_src_vf(i)%sf(l, k, j) = & - flux_src_rsz_vf(j, k, l, i) - end do - end do - end do - end do - - end if - elseif (norm_dir == 1) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, sys_size - do l = is3%beg, is3%end - do k = is2%beg, is2%end - do j = is1%beg, is1%end - flux_vf(i)%sf(j, k, l) = & - flux_rsx_vf(j, k, l, i) + do i = 1, sys_size + flux_vf(i)%sf(l, k, j) = & + flux_rsz_vf(j, k, l, i) + if (grid_geometry == 3) then + flux_gsrc_vf(i)%sf(l, k, j) = & + flux_gsrc_rsz_vf(j, k, l, i) + end if end do end do end do end do + elseif (norm_dir == 1) then !$acc parallel loop collapse(3) gang vector default(present) do l = is3%beg, is3%end do k = is2%beg, is2%end do j = is1%beg, is1%end flux_src_vf(advxb)%sf(j, k, l) = & flux_src_rsx_vf(j, k, l, advxb) + do i = 1, sys_size + flux_vf(i)%sf(j, k, l) = & + flux_rsx_vf(j, k, l, i) + end do end do end do end do + end if - if (riemann_solver == 1 .or. riemann_solver == 4) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = advxb + 1, advxe - do l = is3%beg, is3%end + if (riemann_solver == 1 .or. riemann_solver == 4) then + !$acc parallel loop collapse(4) gang vector default(present) + do i = advxb + 1, advxe + do l = is3%beg, is3%end + do j = is1%beg, is1%end do k = is2%beg, is2%end - do j = is1%beg, is1%end + if (norm_dir == 2) then + flux_src_vf(i)%sf(k, j, l) = & + flux_src_rsy_vf(j, k, l, i) + else if (norm_dir == 3) then + flux_src_vf(i)%sf(l, k, j) = & + flux_src_rsz_vf(j, k, l, i) + else if (norm_dir == 1) then flux_src_vf(i)%sf(j, k, l) = & flux_src_rsx_vf(j, k, l, i) - end do + end if end do end do end do - end if + end do end if end subroutine s_finalize_riemann_solver From 3bfa9f455f0c3a9dd51f0b64a516958ef4b49433 Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Sun, 29 Jun 2025 03:25:19 -0400 Subject: [PATCH 3/9] fixed few errors --- src/simulation/m_riemann_solvers.fpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 8442dad6e..90bbaf03e 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -329,7 +329,7 @@ contains real(wp) :: alpha_L_sum, alpha_R_sum real(wp) :: zcoef, pcorr !< low Mach number correction - type(riemann_states) :: c_fast, pres_mag + type(riemann_states) :: c_fast, pres_mag, vel type(riemann_states_vec3) :: B type(riemann_states) :: Ga ! Gamma (Lorentz factor) @@ -1209,6 +1209,8 @@ contains integer :: i, j, k, l, q !< Generic loop iterators integer :: idx1, idxi + type(riemann_states) :: c_fast, vel + integer :: loop_end ! Populating the buffers of the left and right Riemann problem ! states variables, based on the choice of boundary conditions From ddaa1b6372f87385608bfb5f8b29f14abe55fe95 Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Sun, 29 Jun 2025 03:32:28 -0400 Subject: [PATCH 4/9] quick fix --- src/simulation/m_riemann_solvers.fpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 90bbaf03e..e315db5d5 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -3279,7 +3279,6 @@ contains target, intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, & dqL_prim_dy_vf, dqR_prim_dy_vf, & dqL_prim_dz_vf, dqR_prim_dz_vf, & - qL_prim_vf, qR_prim_vf type(scalar_field), & dimension(:), & pointer :: dqL_prim_d_vf, dqR_prim_d_vf From 12eabf4eced5a95cf7d29e2e3c7121c12d127df5 Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Sun, 29 Jun 2025 18:28:32 -0400 Subject: [PATCH 5/9] simple fix --- src/simulation/m_riemann_solvers.fpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index e315db5d5..527ab4818 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -3279,9 +3279,8 @@ contains target, intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, & dqL_prim_dy_vf, dqR_prim_dy_vf, & dqL_prim_dz_vf, dqR_prim_dz_vf, & - type(scalar_field), & - dimension(:), & - pointer :: dqL_prim_d_vf, dqR_prim_d_vf + + type(scalar_field), allocatable, dimension(:), pointer:: dqL_prim_d_vf, dqR_prim_d_vf integer :: end_val, bc_beg, bc_end From 4e2a439f684ae573d5d196f34b7915b64a438124 Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Sun, 29 Jun 2025 18:43:21 -0400 Subject: [PATCH 6/9] removed & --- src/simulation/m_riemann_solvers.fpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 527ab4818..20ed7ea20 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -161,20 +161,20 @@ contains flux_gsrc_vf, & norm_dir, ix, iy, iz) - real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(INOUT) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf + real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - type(scalar_field), allocatable, dimension(:), intent(INOUT) :: qL_prim_vf, qR_prim_vf + type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf type(scalar_field), & allocatable, dimension(:), & - intent(INOUT) :: dqL_prim_dx_vf, dqR_prim_dx_vf, & + intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, & dqL_prim_dy_vf, dqR_prim_dy_vf, & dqL_prim_dz_vf, dqR_prim_dz_vf type(scalar_field), & dimension(sys_size), & - intent(INOUT) :: flux_vf, flux_src_vf, flux_gsrc_vf + intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf integer, intent(IN) :: norm_dir @@ -225,7 +225,7 @@ contains type(scalar_field), & dimension(sys_size), & - intent(INOUT) :: flux_src_vf + intent(inout) :: flux_src_vf integer, intent(IN) :: norm_dir @@ -3278,7 +3278,7 @@ contains allocatable, dimension(:), & target, intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, & dqL_prim_dy_vf, dqR_prim_dy_vf, & - dqL_prim_dz_vf, dqR_prim_dz_vf, & + dqL_prim_dz_vf, dqR_prim_dz_vf type(scalar_field), allocatable, dimension(:), pointer:: dqL_prim_d_vf, dqR_prim_d_vf From b684a4ea389b74a8d4483f0a29be65bbe638cad0 Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Sun, 29 Jun 2025 18:47:26 -0400 Subject: [PATCH 7/9] remove allocatable from pointers --- src/simulation/m_riemann_solvers.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 20ed7ea20..82cb39bbd 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -3280,7 +3280,7 @@ contains dqL_prim_dy_vf, dqR_prim_dy_vf, & dqL_prim_dz_vf, dqR_prim_dz_vf - type(scalar_field), allocatable, dimension(:), pointer:: dqL_prim_d_vf, dqR_prim_d_vf + type(scalar_field), dimension(:), pointer:: dqL_prim_d_vf, dqR_prim_d_vf integer :: end_val, bc_beg, bc_end From 6fb8c65bfd0991995a1163192dbcd840643a426c Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Sun, 29 Jun 2025 19:15:46 -0400 Subject: [PATCH 8/9] m_riemann_solvers cleanup & format/lint --- src/simulation/m_riemann_solvers.fpp | 20 +------------------- 1 file changed, 1 insertion(+), 19 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 82cb39bbd..dc95953ac 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -687,8 +687,6 @@ contains tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) - xi_M = (5.e-1_wp + sign(5.e-1_wp, s_L)) & + (5.e-1_wp - sign(5.e-1_wp, s_L)) & *(5.e-1_wp + sign(5.e-1_wp, s_R)) @@ -1445,10 +1443,6 @@ contains tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) - ! follows Einfeldt et al. - ! s_M/P = min/max(0.,s_L/R) - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) - ! goes with q_star_L/R = xi_L/R * (variable) ! xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) xi_L = (s_L - vel_L(idx1))/(s_L - s_S) @@ -1718,10 +1712,6 @@ contains tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) - ! follows Einfeldt et al. - ! s_M/P = min/max(0.,s_L/R) - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) - ! goes with q_star_L/R = xi_L/R * (variable) ! xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) @@ -2120,10 +2110,6 @@ contains tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) - ! follows Einfeldt et al. - ! s_M/P = min/max(0.,s_L/R) - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) - ! goes with q_star_L/R = xi_L/R * (variable) ! xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) @@ -2562,10 +2548,6 @@ contains tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) - ! follows Einfeldt et al. - ! s_M/P = min/max(0.,s_L/R) - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) - ! goes with q_star_L/R = xi_L/R * (variable) ! xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) xi_L = (s_L - vel_L(idx1))/(s_L - s_S) @@ -3280,7 +3262,7 @@ contains dqL_prim_dy_vf, dqR_prim_dy_vf, & dqL_prim_dz_vf, dqR_prim_dz_vf - type(scalar_field), dimension(:), pointer:: dqL_prim_d_vf, dqR_prim_d_vf + type(scalar_field), dimension(:), pointer :: dqL_prim_d_vf, dqR_prim_d_vf integer :: end_val, bc_beg, bc_end From dcfc3e02ab23961ec60c98a4556a2ac04578b99f Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Fri, 11 Jul 2025 17:12:51 -0400 Subject: [PATCH 9/9] lint-fix --- src/common/m_variables_conversion.fpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index c150feaf2..3626c90bc 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -1743,11 +1743,11 @@ contains if (wave_speeds == 1) then if (elasticity) then s_L = min(vel_L(idx) - sqrt(c_L*c_L + & - (((4_wp*G_L)/3_wp) + tau_e_L(idx_tau))/rho_L), vel_R(idx) - sqrt(c_R*c_R + & - (((4_wp*G_R)/3_wp) + tau_e_R(idx_tau))/rho_R)) + (((4._wp*G_L)/3._wp) + tau_e_L(idx_tau))/rho_L), vel_R(idx) - sqrt(c_R*c_R + & + (((4._wp*G_R)/3._wp) + tau_e_R(idx_tau))/rho_R)) s_R = max(vel_R(idx) + sqrt(c_R*c_R + & - (((4_wp*G_R)/3_wp) + tau_e_R(idx_tau))/rho_R), vel_L(idx) + sqrt(c_L*c_L + & - (((4_wp*G_L)/3_wp) + tau_e_L(idx_tau))/rho_L)) + (((4._wp*G_R)/3._wp) + tau_e_R(idx_tau))/rho_R), vel_L(idx) + sqrt(c_L*c_L + & + (((4._wp*G_L)/3._wp) + tau_e_L(idx_tau))/rho_L)) s_S = (pres_R - tau_e_R(idx_tau) - pres_L + & tau_e_L(idx_tau) + rho_L*vel_L(idx)*(s_L - vel_L(idx)) - & rho_R*vel_R(idx)*(s_R - vel_R(idx)))/(rho_L*(s_L - vel_L(idx)) - &