Refactor for new variable conversion and modularize sound speed about mfc

[sbryngelson]

For example, see this

MFC/src/simulation/m_data_output.f90

Lines 1645 to 2048 in 83528f1

	subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag) ! -----------
	integer, intent(IN) :: t_step
	type(scalar_field), dimension(sys_size), intent(IN) :: q_cons_vf
	real(kind(0d0)), dimension(0:m, 0:n, 0:p), intent(IN) :: accel_mag
	real(kind(0d0)), dimension(-1:m) :: distx
	real(kind(0d0)), dimension(-1:n) :: disty
	real(kind(0d0)), dimension(-1:p) :: distz
	! The cell-averaged partial densities, density, velocity, pressure,
	! volume fractions, specific heat ratio function, liquid stiffness
	! function, and sound speed.
	real(kind(0d0)) :: lit_gamma, nbub
	real(kind(0d0)) :: rho
	real(kind(0d0)), dimension(num_dims) :: vel
	real(kind(0d0)) :: pres
	real(kind(0d0)) :: ptilde
	real(kind(0d0)) :: ptot
	real(kind(0d0)) :: alf
	real(kind(0d0)) :: alfgr
	real(kind(0d0)), dimension(num_fluids) :: alpha
	real(kind(0d0)) :: gamma
	real(kind(0d0)) :: pi_inf
	real(kind(0d0)) :: c
	real(kind(0d0)) :: M00, M10, M01, M20, M11, M02
	real(kind(0d0)) :: varR, varV
	real(kind(0d0)), dimension(Nb) :: nR, R, nRdot, Rdot
	real(kind(0d0)) :: nR3
	real(kind(0d0)) :: accel
	real(kind(0d0)) :: int_pres
	real(kind(0d0)) :: max_pres
	real(kind(0d0)), dimension(2) :: Re
	real(kind(0d0)) :: E_e
	real(kind(0d0)), dimension(6) :: tau_e
	integer :: i, j, k, l, s !< Generic loop iterator
	real(kind(0d0)) :: nondim_time !< Non-dimensional time
	real(kind(0d0)) :: tmp !<
	!! Temporary variable to store quantity for mpi_allreduce
	real(kind(0d0)) :: blkmod1, blkmod2 !<
	!! Fluid bulk modulus for Woods mixture sound speed
	integer :: npts !< Number of included integral points
	real(kind(0d0)) :: rad, thickness !< For integral quantities
	logical :: trigger !< For integral quantities
	! Non-dimensional time calculation
	if (time_stepper == 23) then
	nondim_time = mytime
	else
	if (t_step_old /= dflt_int) then
	nondim_time = real(t_step + t_step_old, kind(0d0))*dt
	else
	nondim_time = real(t_step, kind(0d0))*dt !*1.d-5/10.0761131451d0
	end if
	end if
	do i = 1, num_probes
	! Zeroing out flow variables for all processors
	rho = 0d0
	do s = 1, num_dims
	vel(s) = 0d0
	end do
	pres = 0d0
	gamma = 0d0
	pi_inf = 0d0
	c = 0d0
	accel = 0d0
	nR = 0d0; R = 0d0
	nRdot = 0d0; Rdot = 0d0
	nbub = 0d0
	M00 = 0d0
	M10 = 0d0
	M01 = 0d0
	M20 = 0d0
	M11 = 0d0
	M02 = 0d0
	varR = 0d0; varV = 0d0
	alf = 0d0
	do s = 1, (num_dims*(num_dims + 1))/2
	tau_e(s) = 0d0
	end do
	! Find probe location in terms of indices on a
	! specific processor
	if (n == 0) then ! 1D simulation
	if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
	do s = -1, m
	distx(s) = x_cb(s) - probe(i)%x
	if (distx(s) < 0d0) distx(s) = 1000d0
	end do
	j = minloc(distx, 1)
	if (j == 1) j = 2 ! Pick first point if probe is at edge
	k = 0
	l = 0
	! Computing/Sharing necessary state variables
	call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, &
	rho, gamma, pi_inf, &
	Re, G, fluid_pp(:)%G)
	do s = 1, num_dims
	vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k, l)/rho
	end do
	if (model_eqns == 4) then
	lit_gamma = 1d0/fluid_pp(1)%gamma + 1d0
	!Tait pressure from density
	pres = (pref + pi_inf)*( &
	(q_cons_vf(1)%sf(j - 2, k, l)/ &
	(rhoref*(1.d0 - q_cons_vf(4)%sf(j - 2, k, l))) &
	)**lit_gamma) &
	- pi_inf
	else if (hypoelasticity) then
	! calculate elastic contribution to Energy
	E_e = 0d0
	do s = stress_idx%beg, stress_idx%end
	if (G > 0) then
	E_e = E_e + ((q_cons_vf(stress_idx%beg)%sf(j - 2, k, l)/rho)**2d0) &
	/(4d0*G)
	! Additional terms in 2D and 3D
	if ((s == stress_idx%beg + 1) .or. &
	(s == stress_idx%beg + 3) .or. &
	(s == stress_idx%beg + 4)) then
	E_e = E_e + ((q_cons_vf(stress_idx%beg)%sf(j - 2, k, l)/rho)**2d0) &
	/(4d0*G)
	end if
	end if
	end do
	tau_e(1) = q_cons_vf(s)%sf(j - 2, k, l)/rho
	pres = ( &
	q_cons_vf(E_idx)%sf(j - 2, k, l) - &
	0.5d0*(q_cons_vf(mom_idx%beg)%sf(j - 2, k, l)**2.d0)/rho - &
	pi_inf - E_e &
	)/gamma
	else if (model_eqns == 2 .and. (bubbles .neqv. .true.)) then
	!Stiffened gas pressure from energy
	pres = ( &
	q_cons_vf(E_idx)%sf(j - 2, k, l) - &
	0.5d0*(q_cons_vf(2)%sf(j - 2, k, l)**2.d0)/q_cons_vf(1)%sf(j - 2, k, l) - &
	pi_inf &
	)/gamma
	else
	!Stiffened gas pressure from energy with bubbles
	pres = ( &
	(q_cons_vf(E_idx)%sf(j - 2, k, l) - &
	0.5d0*(q_cons_vf(mom_idx%beg)%sf(j - 2, k, l)**2.d0)/rho)/ &
	(1.d0 - q_cons_vf(alf_idx)%sf(j - 2, k, l)) - &
	pi_inf &
	)/gamma
	end if
	if (bubbles) then
	alf = q_cons_vf(alf_idx)%sf(j - 2, k, l)
	if (num_fluids == 3) then
	alfgr = q_cons_vf(alf_idx - 1)%sf(j - 2, k, l)
	end if
	do s = 1, nb
	nR(s) = q_cons_vf(bub_idx%rs(s))%sf(j - 2, k, l)
	nRdot(s) = q_cons_vf(bub_idx%vs(s))%sf(j - 2, k, l)
	end do
	!call s_comp_n_from_cons(alf, nR, nbub)
	nR3 = 0d0
	do s = 1, nb
	nR3 = nR3 + weight(s)*(nR(s)**3d0)
	end do
	nbub = DSQRT((4.d0*pi/3.d0)*nR3/alf)
	if (DEBUG) print *, 'In probe, nbub: ', nbub
	if (qbmm) then
	M00 = q_cons_vf(bub_idx%moms(1, 1))%sf(j - 2, k, l)/nbub
	M10 = q_cons_vf(bub_idx%moms(1, 2))%sf(j - 2, k, l)/nbub
	M01 = q_cons_vf(bub_idx%moms(1, 3))%sf(j - 2, k, l)/nbub
	M20 = q_cons_vf(bub_idx%moms(1, 4))%sf(j - 2, k, l)/nbub
	M11 = q_cons_vf(bub_idx%moms(1, 5))%sf(j - 2, k, l)/nbub
	M02 = q_cons_vf(bub_idx%moms(1, 6))%sf(j - 2, k, l)/nbub
	M10 = M10/M00
	M01 = M01/M00
	M20 = M20/M00
	M11 = M11/M00
	M02 = M02/M00
	varR = M20 - M10**2d0
	varV = M02 - M01**2d0
	end if
	R(:) = nR(:)/nbub
	Rdot(:) = nRdot(:)/nbub
	ptilde = ptil(j - 2, k, l)
	ptot = pres - ptilde
	end if
	! Compute mixture sound speed
	if (alt_soundspeed) then
	do s = 1, num_fluids
	alpha(s) = q_cons_vf(E_idx + s)%sf(j - 2, k, l)
	end do
	blkmod1 = ((fluid_pp(1)%gamma + 1d0)*pres + &
	fluid_pp(1)%pi_inf)/fluid_pp(1)%gamma
	blkmod2 = ((fluid_pp(2)%gamma + 1d0)*pres + &
	fluid_pp(2)%pi_inf)/fluid_pp(2)%gamma
	c = (1d0/(rho*(alpha(1)/blkmod1 + alpha(2)/blkmod2)))
	else
	c = (((gamma + 1d0)*pres + pi_inf)/(gamma*rho))
	end if
	if (mixture_err .and. c < 0d0) then
	c = sgm_eps
	else
	c = sqrt(c)
	end if
	accel = accel_mag(j - 2, k, l)
	end if
	elseif (p == 0) then ! 2D simulation
	if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
	if ((probe(i)%y >= y_cb(-1)) .and. (probe(i)%y <= y_cb(n))) then
	do s = -1, m
	distx(s) = x_cb(s) - probe(i)%x
	if (distx(s) < 0d0) distx(s) = 1000d0
	end do
	do s = -1, n
	disty(s) = y_cb(s) - probe(i)%y
	if (disty(s) < 0d0) disty(s) = 1000d0
	end do
	j = minloc(distx, 1)
	k = minloc(disty, 1)
	if (j == 1) j = 2 ! Pick first point if probe is at edge
	if (k == 1) k = 2 ! Pick first point if probe is at edge
	l = 0
	! Computing/Sharing necessary state variables
	call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l, &
	rho, gamma, pi_inf, &
	Re, G, fluid_pp(:)%G)
	do s = 1, num_dims
	vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k - 2, l)/rho
	end do
	if (model_eqns == 4) then
	lit_gamma = 1d0/fluid_pp(1)%gamma + 1d0
	!Tait pressure from density
	pres = (pref + pi_inf)*( &
	(q_cons_vf(1)%sf(j - 2, k - 2, l)/ &
	(rhoref*(1.d0 - q_cons_vf(4)%sf(j - 2, k - 2, l))) &
	)**lit_gamma) &
	- pi_inf
	else if (hypoelasticity) then
	! calculate elastic contribution to Energy
	E_e = 0d0
	do s = stress_idx%beg, stress_idx%end
	if (G > 0) then
	E_e = E_e + ((q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l)/rho)**2d0) &
	/(4d0*G)
	! Additional terms in 2D and 3D
	if ((s == stress_idx%beg + 1) .or. &
	(s == stress_idx%beg + 3) .or. &
	(s == stress_idx%beg + 4)) then
	E_e = E_e + ((q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l)/rho)**2d0) &
	/(4d0*G)
	end if
	end if
	end do
	do s = 1, 3
	tau_e(s) = q_cons_vf(s)%sf(j - 2, k - 2, l)/rho
	end do
	pres = ( &
	q_cons_vf(E_idx)%sf(j - 2, k - 2, l) - &
	0.5d0*(q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l)**2.d0)/rho - &
	pi_inf - E_e &
	)/gamma
	else if (model_eqns == 2 .and. (bubbles .neqv. .true.)) then
	!Stiffened gas pressure from energy
	pres = ( &
	q_cons_vf(E_idx)%sf(j - 2, k - 2, l) - &
	0.5d0*((q_cons_vf(2)%sf(j - 2, k - 2, l)**2.d0 + &
	q_cons_vf(3)%sf(j - 2, k - 2, l)**2.d0)/q_cons_vf(1)%sf(j - 2, k - 2, l)) - &
	pi_inf &
	)/gamma
	else
	!Stiffened gas pressure from energy with bubbles
	pres = ( &
	(q_cons_vf(E_idx)%sf(j - 2, k - 2, l) - &
	0.5d0*(q_cons_vf(2)%sf(j - 2, k - 2, l)**2.d0 + q_cons_vf(3)%sf(j - 2, k - 2, l)**2.d0) &
	/q_cons_vf(1)%sf(j - 2, k - 2, l))/ &
	(1.d0 - q_cons_vf(alf_idx)%sf(j - 2, k - 2, l)) - &
	pi_inf &
	)/gamma
	end if
	if (bubbles) then
	alf = q_cons_vf(alf_idx)%sf(j - 2, k - 2, l)
	do s = 1, nb
	nR(s) = q_cons_vf(bub_idx%rs(s))%sf(j - 2, k - 2, l)
	nRdot(s) = q_cons_vf(bub_idx%vs(s))%sf(j - 2, k - 2, l)
	end do
	!call s_comp_n_from_cons(alf, nR, nbub)
	nR3 = 0d0
	do s = 1, nb
	nR3 = nR3 + weight(s)*(nR(s)**3d0)
	end do
	nbub = DSQRT((4.d0*pi/3.d0)*nR3/alf)
	R(:) = nR(:)/nbub
	Rdot(:) = nRdot(:)/nbub
	end if
	! Compute mixture sound speed
	if (alt_soundspeed) then
	do s = 1, num_fluids
	alpha(s) = q_cons_vf(E_idx + s)%sf(j - 2, k - 2, l)
	end do
	blkmod1 = ((fluid_pp(1)%gamma + 1d0)*pres + &
	fluid_pp(1)%pi_inf)/fluid_pp(1)%gamma
	blkmod2 = ((fluid_pp(2)%gamma + 1d0)*pres + &
	fluid_pp(2)%pi_inf)/fluid_pp(2)%gamma
	c = (1d0/(rho*(alpha(1)/blkmod1 + alpha(2)/blkmod2)))
	else
	c = (((gamma + 1d0)*pres + pi_inf)/(gamma*rho))
	end if
	if (mixture_err .and. c < 0d0) then
	c = sgm_eps
	else
	c = sqrt(c)
	end if
	accel = accel_mag(j - 2, k - 2, l)
	end if
	end if
	else ! 3D simulation
	if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
	if ((probe(i)%y >= y_cb(-1)) .and. (probe(i)%y <= y_cb(n))) then
	if ((probe(i)%z >= z_cb(-1)) .and. (probe(i)%z <= z_cb(p))) then
	do s = -1, m
	distx(s) = x_cb(s) - probe(i)%x
	if (distx(s) < 0d0) distx(s) = 1000d0
	end do
	do s = -1, n
	disty(s) = y_cb(s) - probe(i)%y
	if (disty(s) < 0d0) disty(s) = 1000d0
	end do
	do s = -1, p
	distz(s) = z_cb(s) - probe(i)%z
	if (distz(s) < 0d0) distz(s) = 1000d0
	end do
	j = minloc(distx, 1)
	k = minloc(disty, 1)
	l = minloc(distz, 1)
	if (j == 1) j = 2 ! Pick first point if probe is at edge
	if (k == 1) k = 2 ! Pick first point if probe is at edge
	if (l == 1) l = 2 ! Pick first point if probe is at edge
	! Computing/Sharing necessary state variables
	call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l - 2, &
	rho, gamma, pi_inf, &
	Re, G, fluid_pp(:)%G)
	do s = 1, num_dims
	vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k - 2, l - 2)/rho
	end do
	pres = (q_cons_vf(E_idx)%sf(j - 2, k - 2, l - 2) - 0.5d0*rho*dot_product(vel, vel) - pi_inf)/gamma
	! Compute mixture sound speed
	if (alt_soundspeed) then
	do s = 1, num_fluids
	alpha(s) = q_cons_vf(E_idx + s)%sf(j - 2, k - 2, l - 2)
	end do
	blkmod1 = ((fluid_pp(1)%gamma + 1d0)*pres + &
	fluid_pp(1)%pi_inf)/fluid_pp(1)%gamma
	blkmod2 = ((fluid_pp(2)%gamma + 1d0)*pres + &
	fluid_pp(2)%pi_inf)/fluid_pp(2)%gamma
	c = (1d0/(rho*(alpha(1)/blkmod1 + alpha(2)/blkmod2)))
	else
	c = (((gamma + 1d0)*pres + pi_inf)/(gamma*rho))
	end if
	if (mixture_err .and. c < 0d0) then
	c = sgm_eps
	else
	c = sqrt(c)
	end if
	accel = accel_mag(j - 2, k - 2, l - 2)
	end if
	end if
	end if
	end if

[sbryngelson]

and this

MFC/src/simulation/m_data_output.f90

Lines 693 to 756 in 83528f1

	do j = 0, m
	call s_convert_to_mixture_variables(q_cons_vf, j, 0, 0, rho, gamma, pi_inf, Re, &
	G, fluid_pp(:)%G)
	lit_gamma = 1d0/gamma + 1d0
	if (((i >= cont_idx%beg) .and. (i <= cont_idx%end)) &
	.or. &
	((i >= adv_idx%beg) .and. (i <= adv_idx%end)) &
	) then
	write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
	else if (i == mom_idx%beg) then !u
	write (2, FMT) x_cb(j), q_cons_vf(mom_idx%beg)%sf(j, 0, 0)/rho
	else if (i == stress_idx%beg) then !tau_e
	write (2, FMT) x_cb(j), q_cons_vf(stress_idx%beg)%sf(j, 0, 0)/rho
	else if (i == E_idx) then !p
	if (model_eqns == 4) then
	!Tait pressure from density
	write (2, FMT) x_cb(j), &
	(pref + pi_inf)*( &
	(q_cons_vf(1)%sf(j, 0, 0)/ &
	(rhoref*(1.d0 - q_cons_vf(4)%sf(j, 0, 0))) &
	)**lit_gamma) &
	- pi_inf
	else if (hypoelasticity) then
	! elastic contribution to energy
	E_e = 0d0
	do k = stress_idx%beg, stress_idx%end
	if (G > 1000) then
	E_e = E_e + ((q_cons_vf(stress_idx%beg)%sf(j, 0, 0)/rho)**2d0) &
	/(4d0*G)
	! Additional terms in 2D and 3D
	if ((k == stress_idx%beg + 1) .or. &
	(k == stress_idx%beg + 3) .or. &
	(k == stress_idx%beg + 4)) then
	E_e = E_e + ((q_cons_vf(stress_idx%beg)%sf(j, 0, 0)/rho)**2d0) &
	/(4d0*G)
	end if
	end if
	end do
	write (2, FMT) x_cb(j), &
	( &
	q_cons_vf(E_idx)%sf(j, 0, 0) - &
	0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, 0, 0)**2.d0)/rho - &
	pi_inf - E_e &
	)/gamma
	else if (model_eqns == 2 .and. (bubbles .neqv. .true.)) then
	!Stiffened gas pressure from energy
	write (2, FMT) x_cb(j), &
	( &
	q_cons_vf(E_idx)%sf(j, 0, 0) - &
	0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, 0, 0)**2.d0)/rho - &
	pi_inf &
	)/gamma
	else
	!Stiffened gas pressure from energy with bubbles
	write (2, FMT) x_cb(j), &
	( &
	(q_cons_vf(E_idx)%sf(j, 0, 0) - &
	0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, 0, 0)**2.d0)/rho)/ &
	(1.d0 - q_cons_vf(alf_idx)%sf(j, 0, 0)) - &
	pi_inf &
	)/gamma
	end if

[anshgupta1234]

Update: s_compute_speed_of_sound has now been implemented on anshgupta1234/fix_48. Now moving on to refactoring for new variable conversion

[anshgupta1234]

This should perhaps be closed now?

[sbryngelson]

closed with #86