Add Q_M-criterion to post_process

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

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

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

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, &

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
 

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):