From a36f8f273ca181c25be98c2ce202c335c3fbe8a5 Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Mon, 30 Jun 2025 10:35:02 -0400 Subject: [PATCH 1/7] HLLC solver refactor --- src/simulation/m_riemann_solvers.fpp | 405 ++++++++++----------------- 1 file changed, 145 insertions(+), 260 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 9c65e51e2..c0efa4000 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -1332,15 +1332,6 @@ contains qL_prim_rs${XYZ}$_vf(j, k, l, i) = max(0._wp, qL_prim_rs${XYZ}$_vf(j, k, l, i)) qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) = min(max(0._wp, qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)), 1._wp) alpha_L_sum = alpha_L_sum + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) - end do - - !$acc loop seq - do i = 1, num_fluids - qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)/max(alpha_L_sum, sgm_eps) - end do - - !$acc loop seq - do i = 1, num_fluids qR_prim_rs${XYZ}$_vf(j + 1, k, l, i) = max(0._wp, qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)) qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) = min(max(0._wp, qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)), 1._wp) alpha_R_sum = alpha_R_sum + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) @@ -1348,6 +1339,7 @@ contains !$acc loop seq do i = 1, num_fluids + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)/max(alpha_L_sum, sgm_eps) qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)/max(alpha_R_sum, sgm_eps) end do end if @@ -1372,91 +1364,76 @@ contains !$acc loop seq do i = 1, 2 Re_L(i) = dflt_real - + Re_R(i) = dflt_real if (Re_size(i) > 0) Re_L(i) = 0._wp - + if (Re_size(i) > 0) Re_R(i) = 0._wp !$acc loop seq do q = 1, Re_size(i) Re_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + Re_idx(i, q))/Res(i, q) & + Re_L(i) - end do - - Re_L(i) = 1._wp/max(Re_L(i), sgm_eps) - - end do - - !$acc loop seq - do i = 1, 2 - Re_R(i) = dflt_real - - if (Re_size(i) > 0) Re_R(i) = 0._wp - - !$acc loop seq - do q = 1, Re_size(i) Re_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + Re_idx(i, q))/Res(i, q) & + Re_R(i) end do - + Re_L(i) = 1._wp/max(Re_L(i), sgm_eps) Re_R(i) = 1._wp/max(Re_R(i), sgm_eps) end do end if - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L + E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R - ! ENERGY ADJUSTMENTS FOR HYPOELASTIC ENERGY - if (hypoelasticity) then - !$acc loop seq - do i = 1, strxe - strxb + 1 - tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) - tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) - end do - G_L = 0._wp; G_R = 0._wp + ! ENERGY ADJUSTMENTS FOR HYPOELASTIC/HYPERELASTIC ENERGY + if (hypoelasticity .or. hyperelasticity) then + G_L = 0_wp; G_R = 0_wp !$acc loop seq do i = 1, num_fluids G_L = G_L + alpha_L(i)*Gs(i) G_R = G_R + alpha_R(i)*Gs(i) end do - !$acc loop seq - do i = 1, strxe - strxb + 1 - ! Elastic contribution to energy if G large enough - if ((G_L > verysmall) .and. (G_R > verysmall)) then - E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) - E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) - ! Additional terms in 2D and 3D - if ((i == 2) .or. (i == 4) .or. (i == 5)) then - E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) - E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) + if (hypoelasticity) then + !$acc loop seq + do i = 1, strxe - strxb + 1 + tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) + tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) + end do + !$acc loop seq + do i = 1, strxe - strxb + 1 + ! Elastic contribution to energy if G large enough + if ((G_L > verysmall) .and. (G_R > verysmall)) then + E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4_wp*G_L) + E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4_wp*G_R) + ! Additional terms in 2D and 3D + if ((i == 2) .or. (i == 4) .or. (i == 5)) then + E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4_wp*G_L) + E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4_wp*G_R) + end if end if + end do + else if (hyperelasticity) then + !$acc loop seq + do i = 1, num_dims + xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) + xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) + end do + G_L = 0_wp; G_R = 0_wp; + !$acc loop seq + do i = 1, num_fluids + ! Mixture left and right shear modulus + G_L = G_L + alpha_L(i)*Gs(i) + G_R = G_R + alpha_R(i)*Gs(i) + end do + ! Elastic contribution to energy if G large enough + if (G_L > verysmall .and. G_R > verysmall) then + E_L = E_L + G_L*qL_prim_rs${XYZ}$_vf(j, k, l, xiend + 1) + E_R = E_R + G_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, xiend + 1) end if - end do - end if - - ! ENERGY ADJUSTMENTS FOR HYPERELASTIC ENERGY - if (hyperelasticity) then - !$acc loop seq - do i = 1, num_dims - xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) - xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) - end do - G_L = 0._wp; G_R = 0._wp; - !$acc loop seq - do i = 1, num_fluids - ! Mixture left and right shear modulus - G_L = G_L + alpha_L(i)*Gs(i) - G_R = G_R + alpha_R(i)*Gs(i) - end do - ! Elastic contribution to energy if G large enough - if (G_L > verysmall .and. G_R > verysmall) then - E_L = E_L + G_L*qL_prim_rs${XYZ}$_vf(j, k, l, xiend + 1) - E_R = E_R + G_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, xiend + 1) + !$acc loop seq + do i = 1, b_size - 1 + tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) + tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) + end do end if - !$acc loop seq - do i = 1, b_size - 1 - tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) - tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) - end do end if H_L = (E_L + pres_L)/rho_L @@ -1725,6 +1702,13 @@ contains do k = is2%beg, is2%end do j = is1%beg, is1%end + ! Initialize all variables + vel_L_rms = 0._wp; vel_R_rms = 0._wp + rho_L = 0._wp; rho_R = 0._wp; + gamma_L = 0._wp; gamma_R = 0._wp; + pi_inf_L = 0._wp; pi_inf_R = 0._wp; + qv_L = 0._wp; qv_R = 0._wp; + !$acc loop seq do i = 1, contxe alpha_rho_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, i) @@ -1735,11 +1719,6 @@ contains do i = 1, num_dims vel_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, contxe + i) vel_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, contxe + i) - end do - - vel_L_rms = 0._wp; vel_R_rms = 0._wp - !$acc loop seq - do i = 1, num_dims vel_L_rms = vel_L_rms + vel_L(i)**2._wp vel_R_rms = vel_R_rms + vel_R(i)**2._wp end do @@ -1748,37 +1727,20 @@ contains do i = 1, num_fluids alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) - end do - - pres_L = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx) - pres_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx) - - rho_L = 0._wp - gamma_L = 0._wp - pi_inf_L = 0._wp - qv_L = 0._wp - !$acc loop seq - do i = 1, num_fluids rho_L = rho_L + alpha_rho_L(i) - gamma_L = gamma_L + alpha_L(i)*gammas(i) - pi_inf_L = pi_inf_L + alpha_L(i)*pi_infs(i) - qv_L = qv_L + alpha_rho_L(i)*qvs(i) - end do - - rho_R = 0._wp - gamma_R = 0._wp - pi_inf_R = 0._wp - qv_R = 0._wp - !$acc loop seq - do i = 1, num_fluids rho_R = rho_R + alpha_rho_R(i) + gamma_L = gamma_L + alpha_L(i)*gammas(i) gamma_R = gamma_R + alpha_R(i)*gammas(i) + pi_inf_L = pi_inf_L + alpha_L(i)*pi_infs(i) pi_inf_R = pi_inf_R + alpha_R(i)*pi_infs(i) + qv_L = qv_L + alpha_rho_L(i)*qvs(i) qv_R = qv_R + alpha_rho_R(i)*qvs(i) end do + + pres_L = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx) + pres_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx) E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L @@ -1870,18 +1832,13 @@ contains (1._wp - dir_flg(dir_idx(i)))* & vel_R(dir_idx(i))) - vel_R(dir_idx(i)))) + & dir_flg(dir_idx(i))*pres_R) - end do - - if (bubbles_euler) then - ! Put p_tilde in - !$acc loop seq - do i = 1, num_dims + if (bubbles_euler) then + ! Put p_tilde in flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) = & - flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) + & - xi_M*(dir_flg(dir_idx(i))*(-1._wp*ptilde_L)) & - + xi_P*(dir_flg(dir_idx(i))*(-1._wp*ptilde_R)) - end do - end if + flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) - & + dir_flg(dir_idx(i))*(xi_M*ptilde_L + xi_P*ptilde_R) + end if + end do flux_rs${XYZ}$_vf(j, k, l, E_idx) = 0._wp @@ -1976,14 +1933,19 @@ contains do k = is2%beg, is2%end do j = is1%beg, is1%end + ! Initialize all variables + vel_L_rms = 0._wp; vel_R_rms = 0._wp + rho_L = 0._wp; rho_R = 0._wp; + gamma_L = 0._wp; gamma_R = 0._wp; + pi_inf_L = 0._wp; pi_inf_R = 0._wp; + qv_L = 0._wp; qv_R = 0._wp; + !$acc loop seq do i = 1, num_fluids alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) end do - vel_L_rms = 0._wp; vel_R_rms = 0._wp - !$acc loop seq do i = 1, num_dims vel_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, contxe + i) @@ -1995,57 +1957,28 @@ contains pres_L = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx) pres_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx) - rho_L = 0._wp - gamma_L = 0._wp - pi_inf_L = 0._wp - qv_L = 0._wp + ! Retain this in the refactor + loop_end = num_fluids + if (.not. mpp_lim .and. num_fluids > 2) loop_end = num_fluids - 1 ! Retain this in the refactor if (mpp_lim .and. (num_fluids > 2)) then !$acc loop seq - do i = 1, num_fluids - rho_L = rho_L + qL_prim_rs${XYZ}$_vf(j, k, l, i) - gamma_L = gamma_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*gammas(i) - pi_inf_L = pi_inf_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*pi_infs(i) - qv_L = qv_L + qL_prim_rs${XYZ}$_vf(j, k, l, i)*qvs(i) - end do - else if (num_fluids > 2) then - !$acc loop seq - do i = 1, num_fluids - 1 + do i = 1, loop_end rho_L = rho_L + qL_prim_rs${XYZ}$_vf(j, k, l, i) gamma_L = gamma_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*gammas(i) pi_inf_L = pi_inf_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*pi_infs(i) qv_L = qv_L + qL_prim_rs${XYZ}$_vf(j, k, l, i)*qvs(i) - end do - else - rho_L = qL_prim_rs${XYZ}$_vf(j, k, l, 1) - gamma_L = gammas(1) - pi_inf_L = pi_infs(1) - qv_L = qvs(1) - end if - - rho_R = 0._wp - gamma_R = 0._wp - pi_inf_R = 0._wp - qv_R = 0._wp - - if (mpp_lim .and. (num_fluids > 2)) then - !$acc loop seq - do i = 1, num_fluids - rho_R = rho_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i) - gamma_R = gamma_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*gammas(i) - pi_inf_R = pi_inf_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*pi_infs(i) - qv_R = qv_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)*qvs(i) - end do - else if (num_fluids > 2) then - !$acc loop seq - do i = 1, num_fluids - 1 rho_R = rho_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i) gamma_R = gamma_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*gammas(i) pi_inf_R = pi_inf_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*pi_infs(i) qv_R = qv_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)*qvs(i) end do else + rho_L = qL_prim_rs${XYZ}$_vf(j, k, l, 1) + gamma_L = gammas(1) + pi_inf_L = pi_infs(1) + qv_L = qvs(1) rho_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, 1) gamma_R = gammas(1) pi_inf_R = pi_infs(1) @@ -2057,38 +1990,25 @@ contains !$acc loop seq do i = 1, 2 Re_L(i) = dflt_real - - if (Re_size(i) > 0) Re_L(i) = 0._wp - + Re_R(i) = dflt_real + if (Re_size(i) > 0) then + Re_L(i) = 0._wp + Re_R(i) = 0._wp + end if !$acc loop seq do q = 1, Re_size(i) Re_L(i) = (1._wp - qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + Re_idx(i, q)))/Res(i, q) & + Re_L(i) - end do - - Re_L(i) = 1._wp/max(Re_L(i), sgm_eps) - - end do - - !$acc loop seq - do i = 1, 2 - Re_R(i) = dflt_real - - if (Re_size(i) > 0) Re_R(i) = 0._wp - - !$acc loop seq - do q = 1, Re_size(i) Re_R(i) = (1._wp - qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + Re_idx(i, q)))/Res(i, q) & + Re_R(i) end do - + Re_L(i) = 1._wp/max(Re_L(i), sgm_eps) Re_R(i) = 1._wp/max(Re_R(i), sgm_eps) end do end if end if E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms H_L = (E_L + pres_L)/rho_L @@ -2448,15 +2368,17 @@ contains !idx1 = 1; if (dir_idx(1) == 2) idx1 = 2; if (dir_idx(1) == 3) idx1 = 3 + vel_L_rms = 0._wp; vel_R_rms = 0._wp + rho_L = 0._wp; rho_R = 0._wp + gamma_L = 0._wp; gamma_R = 0._wp + pi_inf_L = 0._wp; pi_inf_R = 0._wp + qv_L = 0._wp; qv_R = 0._wp + alpha_L_sum = 0._wp; alpha_R_sum = 0._wp + !$acc loop seq do i = 1, num_fluids alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) - end do - - vel_L_rms = 0._wp; vel_R_rms = 0._wp - !$acc loop seq - do i = 1, num_dims vel_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, contxe + i) vel_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, contxe + i) vel_L_rms = vel_L_rms + vel_L(i)**2._wp @@ -2466,19 +2388,6 @@ contains pres_L = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx) pres_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx) - rho_L = 0._wp - gamma_L = 0._wp - pi_inf_L = 0._wp - qv_L = 0._wp - - rho_R = 0._wp - gamma_R = 0._wp - pi_inf_R = 0._wp - qv_R = 0._wp - - alpha_L_sum = 0._wp - alpha_R_sum = 0._wp - ! Change this by splitting it into the cases ! present in the bubbles_euler if (mpp_lim) then @@ -2487,22 +2396,13 @@ contains qL_prim_rs${XYZ}$_vf(j, k, l, i) = max(0._wp, qL_prim_rs${XYZ}$_vf(j, k, l, i)) qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) = min(max(0._wp, qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)), 1._wp) alpha_L_sum = alpha_L_sum + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) - end do - - !$acc loop seq - do i = 1, num_fluids - qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)/max(alpha_L_sum, sgm_eps) - end do - - !$acc loop seq - do i = 1, num_fluids qR_prim_rs${XYZ}$_vf(j + 1, k, l, i) = max(0._wp, qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)) qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) = min(max(0._wp, qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)), 1._wp) alpha_R_sum = alpha_R_sum + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) end do - !$acc loop seq do i = 1, num_fluids + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)/max(alpha_L_sum, sgm_eps) qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)/max(alpha_R_sum, sgm_eps) end do end if @@ -2524,31 +2424,19 @@ contains !$acc loop seq do i = 1, 2 Re_L(i) = dflt_real - - if (Re_size(i) > 0) Re_L(i) = 0._wp - + Re_R(i) = dflt_real + if (Re_size(i) > 0) then + Re_L(i) = 0._wp + Re_R(i) = 0._wp + end if !$acc loop seq do q = 1, Re_size(i) Re_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + Re_idx(i, q))/Res(i, q) & + Re_L(i) - end do - - Re_L(i) = 1._wp/max(Re_L(i), sgm_eps) - - end do - - !$acc loop seq - do i = 1, 2 - Re_R(i) = dflt_real - - if (Re_size(i) > 0) Re_R(i) = 0._wp - - !$acc loop seq - do q = 1, Re_size(i) Re_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + Re_idx(i, q))/Res(i, q) & + Re_R(i) end do - + Re_L(i) = 1._wp/max(Re_L(i), sgm_eps) Re_R(i) = 1._wp/max(Re_R(i), sgm_eps) end do end if @@ -2612,60 +2500,57 @@ contains H_R = (E_R + pres_R)/rho_R end if - ! ENERGY ADJUSTMENTS FOR HYPOELASTIC ENERGY - if (hypoelasticity) then - !$acc loop seq - do i = 1, strxe - strxb + 1 - tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) - tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) - end do - G_L = 0._wp - G_R = 0._wp + ! ENERGY ADJUSTMENTS FOR HYPOELASTIC/HYPERELASTIC ENERGY + if (hypoelasticity .or. hyperelasticity) then + G_L = 0_wp; G_R = 0_wp !$acc loop seq do i = 1, num_fluids G_L = G_L + alpha_L(i)*Gs(i) G_R = G_R + alpha_R(i)*Gs(i) end do - !$acc loop seq - do i = 1, strxe - strxb + 1 - ! Elastic contribution to energy if G large enough - if ((G_L > verysmall) .and. (G_R > verysmall)) then - E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) - E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) - ! Additional terms in 2D and 3D - if ((i == 2) .or. (i == 4) .or. (i == 5)) then - E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) - E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) + if (hypoelasticity) then + !$acc loop seq + do i = 1, strxe - strxb + 1 + tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) + tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) + end do + !$acc loop seq + do i = 1, strxe - strxb + 1 + ! Elastic contribution to energy if G large enough + if ((G_L > verysmall) .and. (G_R > verysmall)) then + E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4_wp*G_L) + E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4_wp*G_R) + ! Additional terms in 2D and 3D + if ((i == 2) .or. (i == 4) .or. (i == 5)) then + E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4_wp*G_L) + E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4_wp*G_R) + end if end if + end do + else if (hyperelasticity) then + !$acc loop seq + do i = 1, num_dims + xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) + xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) + end do + G_L = 0_wp; G_R = 0_wp; + !$acc loop seq + do i = 1, num_fluids + ! Mixture left and right shear modulus + G_L = G_L + alpha_L(i)*Gs(i) + G_R = G_R + alpha_R(i)*Gs(i) + end do + ! Elastic contribution to energy if G large enough + if (G_L > verysmall .and. G_R > verysmall) then + E_L = E_L + G_L*qL_prim_rs${XYZ}$_vf(j, k, l, xiend + 1) + E_R = E_R + G_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, xiend + 1) end if - end do - end if - - ! ENERGY ADJUSTMENTS FOR HYPERELASTIC ENERGY - if (hyperelasticity) then - !$acc loop seq - do i = 1, num_dims - xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) - xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) - end do - G_L = 0._wp - G_R = 0._wp - !$acc loop seq - do i = 1, num_fluids - ! Mixture left and right shear modulus - G_L = G_L + alpha_L(i)*Gs(i) - G_R = G_R + alpha_R(i)*Gs(i) - end do - ! Elastic contribution to energy if G large enough - if (G_L > verysmall .and. G_R > verysmall) then - E_L = E_L + G_L*qL_prim_rs${XYZ}$_vf(j, k, l, xiend + 1) - E_R = E_R + G_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, xiend + 1) + !$acc loop seq + do i = 1, b_size - 1 + tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) + tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) + end do end if - !$acc loop seq - do i = 1, b_size - 1 - tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) - tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) - end do end if H_L = (E_L + pres_L)/rho_L From d5dcadd67ce12fa36de8a5c665c12852a3b468bd Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Mon, 30 Jun 2025 10:40:36 -0400 Subject: [PATCH 2/7] declare loop_end --- src/simulation/m_riemann_solvers.fpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index c0efa4000..2b4cd393f 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -1259,6 +1259,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 05bd3a661291cce4dd273e7bd92e6c2560d05dec Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Mon, 30 Jun 2025 10:50:53 -0400 Subject: [PATCH 3/7] lint/format --- src/simulation/m_riemann_solvers.fpp | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 2b4cd393f..ac006e106 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -1381,7 +1381,6 @@ contains end do end if - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R @@ -1706,11 +1705,10 @@ contains ! Initialize all variables vel_L_rms = 0._wp; vel_R_rms = 0._wp - rho_L = 0._wp; rho_R = 0._wp; - gamma_L = 0._wp; gamma_R = 0._wp; - pi_inf_L = 0._wp; pi_inf_R = 0._wp; - qv_L = 0._wp; qv_R = 0._wp; - + rho_L = 0._wp; rho_R = 0._wp; + gamma_L = 0._wp; gamma_R = 0._wp; + pi_inf_L = 0._wp; pi_inf_R = 0._wp; + qv_L = 0._wp; qv_R = 0._wp; !$acc loop seq do i = 1, contxe alpha_rho_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, i) @@ -1738,7 +1736,7 @@ contains qv_L = qv_L + alpha_rho_L(i)*qvs(i) qv_R = qv_R + alpha_rho_R(i)*qvs(i) end do - + pres_L = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx) pres_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx) @@ -1837,8 +1835,8 @@ contains if (bubbles_euler) then ! Put p_tilde in flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) = & - flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) - & - dir_flg(dir_idx(i))*(xi_M*ptilde_L + xi_P*ptilde_R) + flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) - & + dir_flg(dir_idx(i))*(xi_M*ptilde_L + xi_P*ptilde_R) end if end do @@ -1937,11 +1935,10 @@ contains ! Initialize all variables vel_L_rms = 0._wp; vel_R_rms = 0._wp - rho_L = 0._wp; rho_R = 0._wp; - gamma_L = 0._wp; gamma_R = 0._wp; - pi_inf_L = 0._wp; pi_inf_R = 0._wp; - qv_L = 0._wp; qv_R = 0._wp; - + rho_L = 0._wp; rho_R = 0._wp; + gamma_L = 0._wp; gamma_R = 0._wp; + pi_inf_L = 0._wp; pi_inf_R = 0._wp; + qv_L = 0._wp; qv_R = 0._wp; !$acc loop seq do i = 1, num_fluids alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) From 3f03dc67ceb70d58bd6f11b0ab6d285061288d75 Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Mon, 30 Jun 2025 23:11:06 -0400 Subject: [PATCH 4/7] quick fixes to HLLC --- src/simulation/m_riemann_solvers.fpp | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index ac006e106..7e82a3764 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -1302,7 +1302,13 @@ contains idx1 = dir_idx(1) + ! Initialize variables vel_L_rms = 0._wp; vel_R_rms = 0._wp + rho_L = 0._wp; rho_R = 0._wp + gamma_L = 0._wp; gamma_R = 0._wp + pi_inf_L = 0._wp; pi_inf_R = 0._wp + qv_L = 0._wp; qv_R = 0._wp + alpha_L_sum = 0._wp; alpha_R_sum = 0._wp !$acc loop seq do i = 1, num_dims @@ -1315,19 +1321,6 @@ contains pres_L = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx) pres_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx) - rho_L = 0._wp - gamma_L = 0._wp - pi_inf_L = 0._wp - qv_L = 0._wp - - rho_R = 0._wp - gamma_R = 0._wp - pi_inf_R = 0._wp - qv_R = 0._wp - - alpha_L_sum = 0._wp - alpha_R_sum = 0._wp - if (mpp_lim) then !$acc loop seq do i = 1, num_fluids @@ -2378,6 +2371,10 @@ contains do i = 1, num_fluids alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) + end do + + !$acc loop seq + do i = 1, num_dims vel_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, contxe + i) vel_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, contxe + i) vel_L_rms = vel_L_rms + vel_L(i)**2._wp From 11bd67d1850b46faccd1b696eb5a0943316ba46a Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Mon, 30 Jun 2025 23:16:05 -0400 Subject: [PATCH 5/7] lint fixes --- src/simulation/m_riemann_solvers.fpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 7e82a3764..cced0cc6f 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -1379,7 +1379,7 @@ contains ! ENERGY ADJUSTMENTS FOR HYPOELASTIC/HYPERELASTIC ENERGY if (hypoelasticity .or. hyperelasticity) then - G_L = 0_wp; G_R = 0_wp + G_L = 0._wp; G_R = 0._wp !$acc loop seq do i = 1, num_fluids G_L = G_L + alpha_L(i)*Gs(i) @@ -1395,12 +1395,12 @@ contains do i = 1, strxe - strxb + 1 ! Elastic contribution to energy if G large enough if ((G_L > verysmall) .and. (G_R > verysmall)) then - E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4_wp*G_L) - E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4_wp*G_R) + E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) + E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) ! Additional terms in 2D and 3D if ((i == 2) .or. (i == 4) .or. (i == 5)) then - E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4_wp*G_L) - E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4_wp*G_R) + E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) + E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) end if end if end do @@ -1410,7 +1410,7 @@ contains xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) end do - G_L = 0_wp; G_R = 0_wp; + G_L = 0._wp; G_R = 0._wp; !$acc loop seq do i = 1, num_fluids ! Mixture left and right shear modulus @@ -2498,7 +2498,7 @@ contains ! ENERGY ADJUSTMENTS FOR HYPOELASTIC/HYPERELASTIC ENERGY if (hypoelasticity .or. hyperelasticity) then - G_L = 0_wp; G_R = 0_wp + G_L = 0._wp; G_R = 0._wp !$acc loop seq do i = 1, num_fluids G_L = G_L + alpha_L(i)*Gs(i) @@ -2514,12 +2514,12 @@ contains do i = 1, strxe - strxb + 1 ! Elastic contribution to energy if G large enough if ((G_L > verysmall) .and. (G_R > verysmall)) then - E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4_wp*G_L) - E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4_wp*G_R) + E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) + E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) ! Additional terms in 2D and 3D if ((i == 2) .or. (i == 4) .or. (i == 5)) then - E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4_wp*G_L) - E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4_wp*G_R) + E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) + E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) end if end if end do @@ -2529,7 +2529,7 @@ contains xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) end do - G_L = 0_wp; G_R = 0_wp; + G_L = 0._wp; G_R = 0._wp; !$acc loop seq do i = 1, num_fluids ! Mixture left and right shear modulus From e904f754af68ab766092693326fad7e62c60267b Mon Sep 17 00:00:00 2001 From: "Mohammed S. Al-Mahrouqi" <145478595+Malmahrouqi3@users.noreply.github.com> Date: Mon, 30 Jun 2025 23:35:35 -0400 Subject: [PATCH 6/7] Update src/simulation/m_riemann_solvers.fpp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- 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 d15d6a7af..4a2f4a79c 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -1162,7 +1162,7 @@ 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 From 77da589e6a33b4533604943e0c8c139810a09710 Mon Sep 17 00:00:00 2001 From: malmahrouqi3 Date: Fri, 11 Jul 2025 17:23:50 -0400 Subject: [PATCH 7/7] format/lint fix --- src/simulation/m_riemann_solvers.fpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 5beab8512..f86407f91 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -1162,7 +1162,7 @@ contains integer :: i, j, k, l, q !< Generic loop iterators integer :: idx1, idxi - + integer :: loop_end ! Populating the buffers of the left and right Riemann problem @@ -1313,7 +1313,7 @@ contains xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) end do - G_L = 0._wp; G_R = 0._wp; + G_L = 0._wp; G_R = 0._wp; !$acc loop seq do i = 1, num_fluids ! Mixture left and right shear modulus @@ -2432,7 +2432,7 @@ contains xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) end do - G_L = 0._wp; G_R = 0._wp; + G_L = 0._wp; G_R = 0._wp; !$acc loop seq do i = 1, num_fluids ! Mixture left and right shear modulus