diff --git a/src/post_process/m_derived_variables.f90 b/src/post_process/m_derived_variables.f90 index 0c9966387c..df0feca0b6 100644 --- a/src/post_process/m_derived_variables.f90 +++ b/src/post_process/m_derived_variables.f90 @@ -26,6 +26,7 @@ module m_derived_variables s_derive_sound_speed, & s_derive_flux_limiter, & s_derive_vorticity_component, & + s_derive_qm, & s_derive_numerical_schlieren_function, & s_finalize_derived_variables_module @@ -569,6 +570,92 @@ subroutine s_derive_vorticity_component(i, q_prim_vf, q_sf) ! ---------- end subroutine s_derive_vorticity_component ! -------------------------- + !> This subroutine gets as inputs the primitive variables. From those + !! inputs, it proceeds to calculate the value of the Q_M + !! function, which are subsequently stored in the derived flow + !! quantity storage variable, q_sf. + !! @param q_prim_vf Primitive variables + !! @param q_sf Q_M + subroutine s_derive_qm(q_prim_vf,q_sf) + type(scalar_field), & + dimension(sys_size), & + intent(IN) :: q_prim_vf + + real(kind(0d0)), & + dimension(-offset_x%beg:m + offset_x%end, & + -offset_y%beg:n + offset_y%end, & + -offset_z%beg:p + offset_z%end), & + intent(INOUT) :: q_sf + + real(kind(0d0)), & + dimension(1:3, 1:3) :: q_jacobian_sf, S, S2, O, O2 + + real(kind(0d0)) :: trS, trS2, trO2, Q, IIS + integer :: j, k, l, r, jj, kk !< Generic loop iterators + + do l = -offset_z%beg, p + offset_z%end + do k = -offset_y%beg, n + offset_y%end + do j = -offset_x%beg, m + offset_x%end + + ! Get velocity gradient tensor + q_jacobian_sf(:, :) = 0d0 + + do r = -fd_number, fd_number + do jj = 1, 3 + ! d()/dx + q_jacobian_sf(jj, 1) = & + q_jacobian_sf(jj, 1)+ & + fd_coeff_x(r, j)* & + q_prim_vf(mom_idx%beg+jj-1)%sf(r + j, k, l) + ! d()/dy + q_jacobian_sf(jj, 2) = & + q_jacobian_sf(jj, 2)+ & + fd_coeff_y(r, k)* & + q_prim_vf(mom_idx%beg+jj-1)%sf(j, r + k, l) + ! d()/dz + q_jacobian_sf(jj, 3) = & + q_jacobian_sf(jj, 3)+ & + fd_coeff_z(r, l)* & + q_prim_vf(mom_idx%beg+jj-1)%sf(j, k, r + l) + end do + end do + + ! Decompose J into asymmetric matrix, S, and a skew-symmetric matrix, O + do jj = 1, 3 + do kk = 1, 3 + S(jj, kk) = 0.5D0* & + (q_jacobian_sf(jj, kk) + q_jacobian_sf(kk, jj)) + O(jj, kk) = 0.5D0* & + (q_jacobian_sf(jj, kk) - q_jacobian_sf(kk, jj)) + end do + end do + + ! Compute S2 = S*S' + do jj = 1, 3 + do kk = 1, 3 + O2(jj, kk) = O(jj,1)*O(kk,1)+ & + O(jj,2)*O(kk,2)+ & + O(jj,3)*O(kk,3) + S2(jj, kk) = S(jj,1)*S(kk,1)+ & + S(jj,2)*S(kk,2)+ & + S(jj,3)*S(kk,3) + end do + end do + + ! Compute Q + Q = 0.5*((O2(1,1)+O2(2,2)+O2(3,3))- & + (S2(1,1)+S2(2,2)+S2(3,3))) + trS = S(1,1)+S(2,2)+S(3,3) + IIS = 0.5*((S(1,1)+S(2,2)+S(3,3))**2- & + (S2(1,1)+S2(2,2)+S2(3,3))) + q_sf(j, k, l) = Q+IIS + + end do + end do + end do + + end subroutine s_derive_qm + !> This subroutine gets as inputs the conservative variables !! and density. From those inputs, it proceeds to calculate !! the values of the numerical Schlieren function, which are diff --git a/src/post_process/m_global_parameters.f90 b/src/post_process/m_global_parameters.f90 index 186b14efe3..ebda428484 100644 --- a/src/post_process/m_global_parameters.f90 +++ b/src/post_process/m_global_parameters.f90 @@ -187,6 +187,7 @@ module m_global_parameters logical :: cons_vars_wrt logical :: c_wrt logical, dimension(3) :: omega_wrt + logical :: qm_wrt logical :: schlieren_wrt !> @} @@ -318,6 +319,7 @@ subroutine s_assign_default_values_to_user_inputs() ! ------------------ cons_vars_wrt = .false. c_wrt = .false. omega_wrt = .false. + qm_wrt = .false. schlieren_wrt = .false. schlieren_alpha = dflt_real @@ -580,7 +582,7 @@ subroutine s_initialize_global_parameters_module() ! ---------------------- buff_size = max(offset_x%beg, offset_x%end, offset_y%beg, & offset_y%end, offset_z%beg, offset_z%end) - if (any(omega_wrt) .or. schlieren_wrt) then + if (any(omega_wrt) .or. schlieren_wrt .or. qm_wrt) then fd_number = max(1, fd_order/2) buff_size = buff_size + fd_number end if diff --git a/src/post_process/m_mpi_proxy.fpp b/src/post_process/m_mpi_proxy.fpp index 980e4e4be2..2a282d9202 100644 --- a/src/post_process/m_mpi_proxy.fpp +++ b/src/post_process/m_mpi_proxy.fpp @@ -165,7 +165,7 @@ contains & 'alt_soundspeed', 'hypoelasticity', 'parallel_io', 'rho_wrt', & & 'coarsen_silo', 'E_wrt', 'pres_wrt', 'gamma_wrt', & & 'heat_ratio_wrt', 'pi_inf_wrt', 'pres_inf_wrt', 'cons_vars_wrt', & - & 'prim_vars_wrt', 'c_wrt', 'schlieren_wrt', 'bubbles', & + & 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt', 'bubbles', & & 'polytropic', 'fourier_decomp', 'polydisperse' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index cfeebdb6c2..53b55bbe15 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -43,7 +43,7 @@ subroutine s_read_input_file() ! --------------------------------------- E_wrt, pres_wrt, alpha_wrt, gamma_wrt, & heat_ratio_wrt, pi_inf_wrt, pres_inf_wrt, & cons_vars_wrt, prim_vars_wrt, c_wrt, & - omega_wrt, schlieren_wrt, schlieren_alpha, & + omega_wrt, qm_wrt, schlieren_wrt, schlieren_alpha, & fd_order, mixture_err, alt_soundspeed, & flux_lim, flux_wrt, cyl_coord, & parallel_io, coarsen_silo, fourier_decomp, & diff --git a/src/post_process/p_main.f90 b/src/post_process/p_main.f90 index 0b29ec199d..7101b4a41e 100644 --- a/src/post_process/p_main.f90 +++ b/src/post_process/p_main.f90 @@ -124,19 +124,19 @@ program p_main call s_write_grid_to_formatted_database_file(t_step) ! Computing centered finite-difference coefficients in x-direction - if (omega_wrt(2) .or. omega_wrt(3) .or. schlieren_wrt) then + if (omega_wrt(2) .or. omega_wrt(3) .or. qm_wrt .or. schlieren_wrt) then call s_compute_finite_difference_coefficients(m, offset_x, x_cc, & fd_coeff_x) end if ! Computing centered finite-difference coefficients in y-direction - if (omega_wrt(1) .or. omega_wrt(3) .or. (n > 0 .and. schlieren_wrt)) then + if (omega_wrt(1) .or. omega_wrt(3) .or. qm_wrt .or. (n > 0 .and. schlieren_wrt)) then call s_compute_finite_difference_coefficients(n, offset_y, y_cc, & fd_coeff_y) end if ! Computing centered finite-difference coefficients in z-direction - if (omega_wrt(1) .or. omega_wrt(2) .or. (p > 0 .and. schlieren_wrt)) then + if (omega_wrt(1) .or. omega_wrt(2) .or. qm_wrt .or. (p > 0 .and. schlieren_wrt)) then call s_compute_finite_difference_coefficients(p, offset_z, z_cc, & fd_coeff_z) end if @@ -420,6 +420,17 @@ program p_main end if ! ---------------------------------------------------------------------- + ! Adding Q_M to the formatted database file ------------------ + if (p > 0 .and. qm_wrt) then + call s_derive_qm(q_prim_vf, q_sf) + + write (varname, '(A)') 'qm' + call s_write_variable_to_formatted_database_file(varname, t_step) + + varname(:) = ' ' + end if + ! ---------------------------------------------------------------------- + ! Adding numerical Schlieren function to formatted database file ------- if (schlieren_wrt) then diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index d2f0857e9b..3be1e0cd42 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -116,7 +116,7 @@ 'mom_wrt', 'vel_wrt', 'flux_lim', 'flux_wrt', 'E_wrt', 'pres_wrt', 'alpha_wrt', 'kappa_wrt', 'gamma_wrt', 'heat_ratio_wrt', 'pi_inf_wrt', 'pres_inf_wrt', 'cons_vars_wrt', 'prim_vars_wrt', 'c_wrt', 'omega_wrt', - 'schlieren_wrt', 'schlieren_alpha', 'fd_order' + 'qm_wrt', 'schlieren_wrt', 'schlieren_alpha', 'fd_order' ] for cmp_id in range(1,3+1):