diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 8b9d0cf04..f86407f91 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -1163,6 +1163,8 @@ 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 ! states variables, based on the choice of boundary conditions @@ -1203,7 +1205,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 @@ -1216,34 +1224,12 @@ 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 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) @@ -1251,6 +1237,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 @@ -1275,91 +1262,75 @@ 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_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 + ! 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 + 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 @@ -1628,6 +1599,12 @@ 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) @@ -1638,11 +1615,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 @@ -1651,37 +1623,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 - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L + 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 @@ -1773,18 +1728,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 @@ -1879,14 +1829,18 @@ 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) @@ -1898,57 +1852,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) @@ -1960,38 +1885,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 @@ -2351,13 +2263,19 @@ 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) @@ -2369,19 +2287,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 @@ -2390,22 +2295,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 @@ -2427,31 +2323,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 @@ -2515,60 +2399,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 + 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