Equation of state is not DRY

[sbryngelson]

The equation of state is used to compute the pressure from the variables. This EOS changes slightly for different models. The EOS implementation is repeated several places in the code, including in m_variables_converison, m_data_output (including serial data output and probe outputs), and likely more. Following the DRY principle (don't repeat yourself), we should fix this. It has also caused multiple false bugs in the past.

[sbryngelson]

This above doesn't even mention the fact that this same problem is repeated in the pre-process and post-process codes. So we effectively do EoS independently in like 8 places (or more).

[sbryngelson]

The probably is probably best addressed (after discussion with @anandrdbz) by creating pointers, assigned in initialize_module that point to specific functions for computing the pressure (associated with the E_idx) from conservative variables as it depends on model_eqns, bubbles, etc. For example, we could just have a pointer called s_compute_pressure that points to, e.g., s_compute_pressure_5eq_bubbles(q_cons_vf,p) for the 5 equation model with sub-grid bubbles (with cons. variables input and pressure output). Then we could have a separate subroutine for each specific case.

[sbryngelson]

This has similarities to new issue MFlowCode/MFC-develop#18

Thus - it also should probably not be addressed until GPU branch is merged.

[anshgupta1234]

I have gone ahead and implemented the s_compute_pressure subroutines and pointers, located in m_initial_condition.fpp. However, because m_initial_condition.fpp uses m_variables_conversion.fpp and m_variables_conversion.fpp uses m_initial_condition.fpp, there is an error being thrown:

~/src/pre_process/m_variables_conversion.f90:19:9:

   19 |     use m_initial_condition
      |         1
Error: 'm_variables_conversion' of module 'm_initial_condition', imported at (1), is also the name of the current program unit

Therefore I am unsure if there's any way I can really implement s_compute_pressure in m_initial_condition. Is there any other module that I can try and implement it in?

[sbryngelson]

I think m_variable_conversion is a more natural home for this kind of thing.

[sbryngelson]

@anshgupta1234 Check out these lines in m_cbc:

MFC/src/simulation/m_cbc.fpp

Lines 786 to 853 in d80602a

	! 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, adv, alpha_rho, 0, k, r)
	else
	call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, adv, alpha_rho, Re_cbc, 0, k, r)
	end if
	E = gamma*pres + pi_inf + 5d-1*rho*vel_K_sum
	H = (E + pres)/rho
	!$acc loop seq
	do i = 1, contxe
	mf(i) = alpha_rho(i)/rho
	end do
	! Compute mixture sound speed
	if (alt_soundspeed) then
	blkmod1 = ((gammas(1) + 1d0)*pres + &
	pi_infs(1))/gammas(1)
	blkmod2 = ((gammas(2) + 1d0)*pres + &
	pi_infs(2))/gammas(2)
	c = (1d0/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
	elseif (model_eqns == 3) then
	c = 0d0
	!$acc loop seq
	do i = 1, num_fluids
	c = c + q_prim_rs${XYZ}$_vf(0, k, r, i + advxb - 1)*(1d0/gammas(i) + 1d0)* &
	(pres + pi_infs(i)/(gammas(i) + 1d0))
	end do
	c = c/rho
	else
	c = ((H - 5d-1*vel_K_sum)/gamma)
	end if
	c = sqrt(c)
	! IF (mixture_err .AND. c < 0d0) THEN
	! c = sgm_eps
	! ELSE
	! c = SQRT(c)
	! END IF
	! ============================================================

I think these can probably be factored out completely. One thing we didn't initially notice is how pervasive the computation of c, the speed of sound is.
This is above, but also throughout the codebase.
We should create a helper function that goes in common/m_var_conversion that computes the speed of sound.

We see something similar in m_riemann_solvers, for example:

MFC/src/simulation/m_riemann_solvers.fpp

Lines 464 to 572 in d80602a

	!$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
	rho_R = 0d0
	gamma_R = 0d0
	pi_inf_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)
	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)
	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
	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

We also see the computation of c coming up again. Each time it takes something like 70 lines:

MFC/src/simulation/m_riemann_solvers.fpp

Lines 639 to 714 in d80602a

	if (mixture_err) then
	if ((H_avg - 5d-1*vel_avg_rms) < 0d0) then
	c_avg = sgm_eps
	else
	c_avg = sqrt((H_avg - 5d-1*vel_avg_rms)/gamma_avg)
	end if
	else
	c_avg = sqrt((H_avg - 5d-1*vel_avg_rms)/gamma_avg)
	end if
	if (alt_soundspeed) then
	blkmod1 = ((gammas(1) + 1d0)*pres_L + &
	pi_infs(1))/gammas(1)
	blkmod2 = ((gammas(2) + 1d0)*pres_L + &
	pi_infs(2))/gammas(2)
	c_L = 1d0/(rho_L*(alpha_L(1)/blkmod1 + alpha_L(2)/blkmod2))
	blkmod1 = ((gammas(1) + 1d0)*pres_R + &
	pi_infs(1))/gammas(1)
	blkmod2 = ((gammas(2) + 1d0)*pres_R + &
	pi_infs(2))/gammas(2)
	c_R = 1d0/(rho_R*(alpha_R(1)/blkmod1 + alpha_R(2)/blkmod2))
	elseif (model_eqns == 3) then
	c_L = 0d0
	c_R = 0d0
	!$acc loop seq
	do i = 1, num_fluids
	c_L = c_L + qL_prim_rs${XYZ}$_vf(j, k, l, i + advxb - 1)*(1d0/gammas(i) + 1d0)* &
	(qL_prim_rs${XYZ}$_vf(j, k, l, E_idx) + pi_infs(i)/(gammas(i) + 1d0))
	c_R = c_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + advxb - 1)*(1d0/gammas(i) + 1d0)* &
	(qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx) + pi_infs(i)/(gammas(i) + 1d0))
	end do
	c_L = c_L/rho_L
	c_R = c_R/rho_R
	elseif ((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles)) then
	! Sound speed for bubble mmixture to order O(\alpha)
	if (mpp_lim .and. (num_fluids > 1)) then
	c_L = (1d0/gamma_L + 1d0)* &
	(pres_L + pi_inf_L)/rho_L
	c_R = (1d0/gamma_R + 1d0)* &
	(pres_R + pi_inf_R)/rho_R
	else
	c_L = &
	(1d0/gamma_L + 1d0)* &
	(pres_L + pi_inf_L)/ &
	(rho_L*(1d0 - alpha_L(num_fluids)))
	c_R = &
	(1d0/gamma_R + 1d0)* &
	(pres_R + pi_inf_R)/ &
	(rho_R*(1d0 - alpha_R(num_fluids)))
	end if
	else
	c_L = ((H_L - 5d-1*vel_L_rms)/gamma_L)
	c_R = ((H_R - 5d-1*vel_R_rms)/gamma_R)
	end if
	if (mixture_err .and. c_L < 0d0) then
	c_L = 100.d0*sgm_eps
	else
	c_L = sqrt(c_L)
	end if
	if (mixture_err .and. c_R < 0d0) then
	c_R = 100.d0*sgm_eps
	else
	c_R = sqrt(c_R)
	end if

[sbryngelson]

I realize this is becoming an omnibus "issue". Let's merge your pull request when you have all of the main variables conversion stuff done, then do speed of sound c after.