Riemann solver enthalpy calculation refactor

[wilfonba]

Once #515 is merged, a new helper function s_compute_enthalpy can be used to reduce the number of lines in m_riemann_solvers by abstracting the enthalpy calculation at the left and right states.

[sbryngelson]

This can presumably also be used anywhere H is computed, which is in several other places (@okBrian this is useful for you as well)

• MFC/src/simulation/m_cbc.fpp

  Lines 769 to 808 in 2651071

  	! Transferring the Primitive Variables =======================
  	!$acc loop seq
  	do i = 1, contxe
  	alpha_rho(i) = q_prim_rs${XYZ}$_vf(0, k, r, i)
  	end do
  	!$acc loop seq
  	do i = 1, num_dims
  	vel(i) = q_prim_rs${XYZ}$_vf(0, k, r, contxe + i)
  	end do
  	vel_K_sum = 0d0
  	!$acc loop seq
  	do i = 1, num_dims
  	vel_K_sum = vel_K_sum + vel(i)**2d0
  	end do
  	pres = q_prim_rs${XYZ}$_vf(0, k, r, E_idx)
  	!$acc loop seq
  	do i = 1, advxe - E_idx
  	adv(i) = q_prim_rs${XYZ}$_vf(0, k, r, E_idx + i)
  	end do
  	if (bubbles) then
  	call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, adv, alpha_rho, Re_cbc, 0, k, r)
  	else
  	call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, adv, alpha_rho, Re_cbc, 0, k, r)
  	end if
  	!$acc loop seq
  	do i = 1, contxe
  	mf(i) = alpha_rho(i)/rho
  	end do
  	E = gamma*pres + pi_inf + 5d-1*rho*vel_K_sum
  	H = (E + pres)/rho

• MFC/src/simulation/m_data_output.fpp

  Lines 252 to 276 in 2651071

  	do i = 1, num_fluids
  	alpha_rho(i) = q_prim_vf(i)%sf(j, k, l)
  	alpha(i) = q_prim_vf(E_idx + i)%sf(j, k, l)
  	end do
  	if (bubbles) then
  	call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l)
  	else
  	call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l)
  	end if
  	do i = 1, num_dims
  	vel(i) = q_prim_vf(contxe + i)%sf(j, k, l)
  	end do
  	vel_sum = 0d0
  	do i = 1, num_dims
  	vel_sum = vel_sum + vel(i)**2d0
  	end do
  	pres = q_prim_vf(E_idx)%sf(j, k, l)
  	E = gamma*pres + pi_inf + 5d-1*rho*vel_sum + qv
  	H = (E + pres)/rho

• MFC/src/simulation/m_riemann_solvers.fpp

  Lines 1270 to 1328 in 2651071

  	!$acc loop seq
  	do i = 1, contxe
  	alpha_rho_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, i)
  	alpha_rho_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, 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)
  	end do
  	vel_L_rms = 0d0; vel_R_rms = 0d0
  	!$acc loop seq
  	do i = 1, num_dims
  	vel_L_rms = vel_L_rms + vel_L(i)**2d0
  	vel_R_rms = vel_R_rms + vel_R(i)**2d0
  	end do
  	!$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
  	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 = 0d0
  	gamma_L = 0d0
  	pi_inf_L = 0d0
  	qv_L = 0d0
  	!$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 = 0d0
  	gamma_R = 0d0
  	pi_inf_R = 0d0
  	qv_R = 0d0
  	!$acc loop seq
  	do i = 1, num_fluids
  	rho_R = rho_R + alpha_rho_R(i)
  	gamma_R = gamma_R + alpha_R(i)*gammas(i)
  	pi_inf_R = pi_inf_R + alpha_R(i)*pi_infs(i)
  	qv_R = qv_R + alpha_rho_R(i)*qvs(i)
  	end do
  	E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + qv_L
  	E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + qv_R
  	H_L = (E_L + pres_L)/rho_L
  	H_R = (E_R + pres_R)/rho_R

• MFC/src/simulation/m_riemann_solvers.fpp

  Lines 1519 to 1636 in 2651071

  	!$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 = 0d0; vel_R_rms = 0d0
  	!$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)**2d0
  	vel_R_rms = vel_R_rms + vel_R(i)**2d0
  	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 = 0d0
  	gamma_L = 0d0
  	pi_inf_L = 0d0
  	qv_L = 0d0
  	! 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
  	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 = 0d0
  	gamma_R = 0d0
  	pi_inf_R = 0d0
  	qv_R = 0d0
  	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_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, 1)
  	gamma_R = gammas(1)
  	pi_inf_R = pi_infs(1)
  	qv_R = qvs(1)
  	end if
  	if (any(Re_size > 0)) then
  	if (num_fluids == 1) then ! Need to consider case with num_fluids >= 2
  	!$acc loop seq
  	do i = 1, 2
  	Re_L(i) = dflt_real
  	if (Re_size(i) > 0) Re_L(i) = 0d0
  	!$acc loop seq
  	do q = 1, Re_size(i)
  	Re_L(i) = (1d0 - 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) = 1d0/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) = 0d0
  	!$acc loop seq
  	do q = 1, Re_size(i)
  	Re_R(i) = (1d0 - qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + Re_idx(i, q)))/Res(i, q) &
  	+ Re_R(i)
  	end do
  	Re_R(i) = 1d0/max(Re_R(i), sgm_eps)
  	end do
  	end if
  	end if
  	E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms
  	E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms
  	H_L = (E_L + pres_L)/rho_L
  	H_R = (E_R + pres_R)/rho_R

• MFC/src/simulation/m_riemann_solvers.fpp

  Lines 361 to 473 in 2651071

  	!$acc loop seq
  	do i = 1, contxe
  	alpha_rho_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, i)
  	alpha_rho_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, 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)
  	end do
  	vel_L_rms = 0d0; vel_R_rms = 0d0
  	!$acc loop seq
  	do i = 1, num_dims
  	vel_L_rms = vel_L_rms + vel_L(i)**2d0
  	vel_R_rms = vel_R_rms + vel_R(i)**2d0
  	end do
  	!$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
  	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 = 0d0
  	gamma_L = 0d0
  	pi_inf_L = 0d0
  	qv_L = 0d0
  	rho_R = 0d0
  	gamma_R = 0d0
  	pi_inf_R = 0d0
  	qv_R = 0d0
  	alpha_L_sum = 0d0
  	alpha_R_sum = 0d0
  	if (mpp_lim) then
  	!$acc loop seq
  	do i = 1, num_fluids
  	alpha_rho_L(i) = max(0d0, alpha_rho_L(i))
  	alpha_L(i) = min(max(0d0, alpha_L(i)), 1d0)
  	alpha_L_sum = alpha_L_sum + alpha_L(i)
  	end do
  	alpha_L = alpha_L/max(alpha_L_sum, sgm_eps)
  	!$acc loop seq
  	do i = 1, num_fluids
  	alpha_rho_R(i) = max(0d0, alpha_rho_R(i))
  	alpha_R(i) = min(max(0d0, alpha_R(i)), 1d0)
  	alpha_R_sum = alpha_R_sum + alpha_R(i)
  	end do
  	alpha_R = alpha_R/max(alpha_R_sum, sgm_eps)
  	end if
  	!$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)
  	rho_R = rho_R + alpha_rho_R(i)
  	gamma_R = gamma_R + alpha_R(i)*gammas(i)
  	pi_inf_R = pi_inf_R + alpha_R(i)*pi_infs(i)
  	qv_R = qv_R + alpha_rho_R(i)*qvs(i)
  	end do
  	if (any(Re_size > 0)) then
  	!$acc loop seq
  	do i = 1, 2
  	Re_L(i) = dflt_real
  	if (Re_size(i) > 0) Re_L(i) = 0d0
  	!$acc loop seq
  	do q = 1, Re_size(i)
  	Re_L(i) = alpha_L(Re_idx(i, q))/Res(i, q) &
  	+ Re_L(i)
  	end do
  	Re_L(i) = 1d0/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) = 0d0
  	!$acc loop seq
  	do q = 1, Re_size(i)
  	Re_R(i) = alpha_R(Re_idx(i, q))/Res(i, q) &
  	+ Re_R(i)
  	end do
  	Re_R(i) = 1d0/max(Re_R(i), sgm_eps)
  	end do
  	end if
  	E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + qv_L
  	E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + qv_R
  	H_L = (E_L + pres_L)/rho_L
  	H_R = (E_R + pres_R)/rho_R

• MFC/src/simulation/m_riemann_solvers.fpp

  Lines 913 to 1028 in 2651071

  	!ME3
  	!$acc parallel loop collapse(3) gang vector default(present) private(vel_L, vel_R, Re_L, Re_R, &
  	!$acc rho_avg, h_avg, gamma_avg, s_L, s_R, s_S, vel_avg_rms, alpha_L, alpha_R)
  	do l = is3%beg, is3%end
  	do k = is2%beg, is2%end
  	do j = is1%beg, is1%end
  	vel_L_rms = 0d0; vel_R_rms = 0d0
  	!$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)**2d0
  	vel_R_rms = vel_R_rms + vel_R(i)**2d0
  	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 = 0d0
  	gamma_L = 0d0
  	pi_inf_L = 0d0
  	qv_L = 0d0
  	rho_R = 0d0
  	gamma_R = 0d0
  	pi_inf_R = 0d0
  	qv_R = 0d0
  	alpha_L_sum = 0d0
  	alpha_R_sum = 0d0
  	if (mpp_lim) then
  	!$acc loop seq
  	do i = 1, num_fluids
  	qL_prim_rs${XYZ}$_vf(j, k, l, i) = max(0d0, qL_prim_rs${XYZ}$_vf(j, k, l, i))
  	qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i) = min(max(0d0, qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)), 1d0)
  	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(0d0, qR_prim_rs${XYZ}$_vf(j + 1, k, l, i))
  	qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i) = min(max(0d0, qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)), 1d0)
  	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
  	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
  	!$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)
  	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)
  	alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, advxb + i - 1)
  	alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, advxb + i - 1)
  	end do
  	if (any(Re_size > 0)) then
  	!$acc loop seq
  	do i = 1, 2
  	Re_L(i) = dflt_real
  	if (Re_size(i) > 0) Re_L(i) = 0d0
  	!$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) = 1d0/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) = 0d0
  	!$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_R(i) = 1d0/max(Re_R(i), sgm_eps)
  	end do
  	end if
  	E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + qv_L
  	E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + qv_R
  	H_L = (E_L + pres_L)/rho_L
  	H_R = (E_R + pres_R)/rho_R

It's also in post-process here (sort of):

MFC/src/post_process/m_start_up.f90

Lines 451 to 452 in 2651071

	H = ((gamma_sf(i, j, k) + 1d0)*pres + &
	pi_inf_sf(i, j, k))/rho_sf(i, j, k)

[sbryngelson]

@ChrisZYJ @wilfonba looking at you guys 👀