Add Liutex to post_process using LAPACK (#970)

Co-authored-by: Hyeoksu Lee <[email protected]>
Co-authored-by: Hyeoksu Lee <[email protected]>
Co-authored-by: Spencer Bryngelson <[email protected]>
Co-authored-by: Hyeoksu Lee <[email protected]>

Author: 5 people

CMakeLists.txt

@@ -391,9 +391,10 @@ HANDLE_SOURCES(syscheck OFF)
 #  * FFTW     (optional) Should be linked with an FFTW-like library (fftw/cufftw),
 #             depending on whether OpenACC is enabled and which compiler is
 #             being used.
+#  * LAPACK   (optional) Should be linked with LAPACK
 
 function(MFC_SETUP_TARGET)
-    cmake_parse_arguments(ARGS "OpenACC;MPI;SILO;HDF5;FFTW" "TARGET" "SOURCES" ${ARGN})
+    cmake_parse_arguments(ARGS "OpenACC;MPI;SILO;HDF5;FFTW;LAPACK" "TARGET" "SOURCES" ${ARGN})
 
     add_executable(${ARGS_TARGET} ${ARGS_SOURCES})
     set(IPO_TARGETS ${ARGS_TARGET})
@@ -461,6 +462,11 @@ function(MFC_SETUP_TARGET)
             endif()
         endif()
 
+				if (ARGS_LAPACK)
+					find_package(LAPACK REQUIRED)
+					target_link_libraries(${a_target} PRIVATE LAPACK::LAPACK)
+				endif()
+
         if (MFC_OpenACC AND ARGS_OpenACC)
             find_package(OpenACC)
 
@@ -547,7 +553,7 @@ endif()
 if (MFC_POST_PROCESS)
     MFC_SETUP_TARGET(TARGET  post_process
                      SOURCES "${post_process_SRCs}"
-                     MPI SILO HDF5 FFTW)
+                     MPI SILO HDF5 FFTW LAPACK)
 
     # -O0 is in response to https://github.com/MFlowCode/MFC-develop/issues/95
     target_compile_options(post_process PRIVATE -O0)

docs/documentation/case.md

@@ -589,6 +589,7 @@ To restart the simulation from $k$-th time step, see [Restarting Cases](running.
 | `omega_wrt(i)`       | Logical | Add the $i$-direction vorticity to the database	 |
 | `schlieren_wrt`      | Logical | Add the numerical schlieren to the database|
 | `qm_wrt`             | Logical | Add the Q-criterion to the database|
+| `liutex_wrt`         | Logical | Add the Liutex to the database|
 | `tau_wrt`            | Logical | Add the elastic stress components to the database|
 | `fd_order`           | Integer | Order of finite differences for computing the vorticity and the numerical Schlieren function [1,2,4] |
 | `schlieren_alpha(i)` | Real    | Intensity of the numerical Schlieren computed via `alpha(i)` |
@@ -628,7 +629,7 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu
 - `output_partial_domain` activates the output of part of the domain specified by `[x,y,z]_output%beg` and `[x,y,z]_output%end`.
 This is useful for large domains where only a portion of the domain is of interest.
 It is not supported when `precision = 1` and `format = 1`. 
-It also cannot be enabled with `flux_wrt`, `heat_ratio_wrt`, `pres_inf_wrt`, `c_wrt`, `omega_wrt`, `ib`, `schlieren_wrt`, or `qm_wrt`.
+It also cannot be enabled with `flux_wrt`, `heat_ratio_wrt`, `pres_inf_wrt`, `c_wrt`, `omega_wrt`, `ib`, `schlieren_wrt`, `qm_wrt`, or 'liutex_wrt'.
 
 ### 8. Acoustic Source {#acoustic-source}
 

docs/documentation/references.md

@@ -69,3 +69,5 @@
 - <a id="Powell94">Powell, K. G. (1994). An approximate Riemann solver for magnetohydrodynamics: (That works in more than one dimension). In Upwind and high-resolution schemes (pp. 570-583). Springer.</a>
 
 - <a id="Cao19">Cao, S., Zhang, Y., Liao, D., Zhong, P., and Wang, K. G. (2019). Shock-induced damage and dynamic fracture in cylindrical bodies submerged in liquid. International Journal of Solids and Structures, 169:55–71. Elsevier.</a>
+
+- <a id="Xu2019">Xu, W., Gao, Y., Deng, Y., Liu, J., and Liu, C. (2019). An explicit expression for the calculation of the Rortex vector. Physics of Fluids, 31(9)..</a>

examples/3D_turb_mixing/README.md

@@ -0,0 +1,4 @@
+# 3D Turbulent Mixing layer (3D)
+
+## Liutex visualization at transitional state
+<img src='result.png' height='MAX_HEIGHT'/>

examples/3D_turb_mixing/case.py

@@ -96,6 +96,7 @@
             "omega_wrt(2)": "T",
             "omega_wrt(3)": "T",
             "qm_wrt": "T",
+            "liutex_wrt": "T",
             # Patch 1
             "patch_icpp(1)%geometry": 9,
             "patch_icpp(1)%x_centroid": Lx / 2.0,

examples/3D_turb_mixing/result.png

@@ -32,6 +32,8 @@ contains
         call s_check_inputs_flux_limiter
         call s_check_inputs_volume_fraction
         call s_check_inputs_vorticity
+        call s_check_inputs_qm
+        call s_check_inputs_liutex
         call s_check_inputs_schlieren
         call s_check_inputs_surface_tension
         call s_check_inputs_no_flow_variables
@@ -49,8 +51,7 @@ contains
     impure subroutine s_check_inputs_partial_domain
         @:PROHIBIT(output_partial_domain .and. format == 1)
         @:PROHIBIT(output_partial_domain .and. precision == 1)
-        @:PROHIBIT(output_partial_domain .and. any([flux_wrt, heat_ratio_wrt, pres_inf_wrt, c_wrt, schlieren_wrt, qm_wrt, ib, any(omega_wrt)]))
-
+        @:PROHIBIT(output_partial_domain .and. any([flux_wrt, heat_ratio_wrt, pres_inf_wrt, c_wrt, schlieren_wrt, qm_wrt, liutex_wrt, ib, any(omega_wrt)]))
         @:PROHIBIT(output_partial_domain .and. (f_is_default(x_output%beg) .or. f_is_default(x_output%end)))
         @:PROHIBIT(output_partial_domain .and. n /= 0 .and. (f_is_default(y_output%beg) .or. f_is_default(y_output%end)))
         @:PROHIBIT(output_partial_domain .and. p /= 0 .and. (f_is_default(z_output%beg) .or. f_is_default(z_output%end)))
@@ -110,6 +111,16 @@ contains
         @:PROHIBIT(any(omega_wrt) .and. fd_order == dflt_int, "fd_order must be set for omega_wrt")
     end subroutine s_check_inputs_vorticity
 
+    !> Checks constraints on Q-criterion parameters
+    impure subroutine s_check_inputs_qm
+        @:PROHIBIT(n == 0 .and. qm_wrt)
+    end subroutine s_check_inputs_qm
+
+    !> Checks constraints on liutex parameters
+    impure subroutine s_check_inputs_liutex
+        @:PROHIBIT(n == 0 .and. liutex_wrt)
+    end subroutine s_check_inputs_liutex
+
     !> Checks constraints on numerical Schlieren parameters
         !! (schlieren_wrt and schlieren_alpha)
     impure subroutine s_check_inputs_schlieren

src/post_process/m_checker.fpp

@@ -29,6 +29,7 @@ module m_derived_variables
  s_derive_flux_limiter, &
  s_derive_vorticity_component, &
  s_derive_qm, &
+ s_derive_liutex, &
  s_derive_numerical_schlieren_function, &
  s_compute_speed_of_sound, &
  s_finalize_derived_variables_module
@@ -80,21 +81,21 @@ contains
         ! s_compute_finite_difference_coefficients.
 
         ! Allocating 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. schlieren_wrt .or. liutex_wrt) then
             allocate (fd_coeff_x(-fd_number:fd_number, &
                                  -offset_x%beg:m + offset_x%end))
         end if
 
         ! Allocating centered finite-difference coefficients in y-direction
-        if (omega_wrt(1) .or. omega_wrt(3) &
+        if (omega_wrt(1) .or. omega_wrt(3) .or. liutex_wrt &
             .or. &
             (n > 0 .and. schlieren_wrt)) then
             allocate (fd_coeff_y(-fd_number:fd_number, &
                                  -offset_y%beg:n + offset_y%end))
         end if
 
         ! Allocating centered finite-difference coefficients in z-direction
-        if (omega_wrt(1) .or. omega_wrt(2) &
+        if (omega_wrt(1) .or. omega_wrt(2) .or. liutex_wrt &
             .or. &
             (p > 0 .and. schlieren_wrt)) then
             allocate (fd_coeff_z(-fd_number:fd_number, &
@@ -556,6 +557,135 @@ contains
 
     end subroutine s_derive_qm
 
+    !> This subroutine gets as inputs the primitive variables. From those
+        !!      inputs, it proceeds to calculate the Liutex vector and its
+        !!      magnitude based on Xu et al. (2019).
+        !!  @param q_prim_vf Primitive variables
+        !!  @param liutex_mag Liutex magnitude
+        !!  @param liutex_axis Liutex axis
+    impure subroutine s_derive_liutex(q_prim_vf, liutex_mag, liutex_axis)
+        integer, parameter :: nm = 3
+        type(scalar_field), &
+            dimension(sys_size), &
+            intent(in) :: q_prim_vf
+
+        real(wp), &
+            dimension(-offset_x%beg:m + offset_x%end, &
+                      -offset_y%beg:n + offset_y%end, &
+                      -offset_z%beg:p + offset_z%end), &
+            intent(out) :: liutex_mag !< Liutex magnitude
+
+        real(wp), &
+            dimension(-offset_x%beg:m + offset_x%end, &
+                      -offset_y%beg:n + offset_y%end, &
+                      -offset_z%beg:p + offset_z%end, nm), &
+            intent(out) :: liutex_axis !< Liutex rigid rotation axis
+
+        character, parameter :: ivl = 'N' !< compute left eigenvectors
+        character, parameter :: ivr = 'V' !< compute right eigenvectors
+        real(wp), dimension(nm, nm) :: vgt !< velocity gradient tensor
+        real(wp), dimension(nm) :: lr, li !< real and imaginary parts of eigenvalues
+        real(wp), dimension(nm, nm) :: vl, vr !< left and right eigenvectors
+        integer, parameter :: lwork = 4*nm !< size of work array (4*nm recommended)
+        real(wp), dimension(lwork) :: work !< work array
+        integer :: info
+
+        real(wp), dimension(nm) :: eigvec !< real eigenvector
+        real(wp) :: eigvec_mag !< magnitude of real eigenvector
+        real(wp) :: omega_proj !< projection of vorticity on real eigenvector
+        real(wp) :: lci !< imaginary part of complex eigenvalue
+        real(wp) :: alpha
+
+        integer :: j, k, l, r, i !< Generic loop iterators
+        integer :: idx
+
+        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 (VGT)
+                    vgt(:, :) = 0._wp
+
+                    do r = -fd_number, fd_number
+                        do i = 1, 3
+                            ! d()/dx
+                            vgt(i, 1) = &
+                                vgt(i, 1) + &
+                                fd_coeff_x(r, j)* &
+                                q_prim_vf(mom_idx%beg + i - 1)%sf(r + j, k, l)
+                            ! d()/dy
+                            vgt(i, 2) = &
+                                vgt(i, 2) + &
+                                fd_coeff_y(r, k)* &
+                                q_prim_vf(mom_idx%beg + i - 1)%sf(j, r + k, l)
+                            ! d()/dz
+                            vgt(i, 3) = &
+                                vgt(i, 3) + &
+                                fd_coeff_z(r, l)* &
+                                q_prim_vf(mom_idx%beg + i - 1)%sf(j, k, r + l)
+                        end do
+                    end do
+
+                    ! Call appropriate LAPACK routine based on precision
+#ifdef MFC_SINGLE_PRECISION
+                    call sgeev(ivl, ivr, nm, vgt, nm, lr, li, vl, nm, vr, nm, work, lwork, info)
+#else
+                    call dgeev(ivl, ivr, nm, vgt, nm, lr, li, vl, nm, vr, nm, work, lwork, info)
+#endif
+
+                    ! Find real eigenvector
+                    idx = 1
+                    do r = 2, 3
+                        if (abs(li(r)) < abs(li(idx))) then
+                            idx = r
+                        end if
+                    end do
+                    eigvec = vr(:, idx)
+
+                    ! Normalize real eigenvector if it is effectively non-zero
+                    eigvec_mag = sqrt(eigvec(1)**2._wp &
+                                      + eigvec(2)**2._wp &
+                                      + eigvec(3)**2._wp)
+                    if (eigvec_mag > sgm_eps) then
+                        eigvec = eigvec/eigvec_mag
+                    else
+                        eigvec = 0._wp
+                    end if
+
+                    ! Compute vorticity projected on the eigenvector
+                    omega_proj = (vgt(3, 2) - vgt(2, 3))*eigvec(1) &
+                                 + (vgt(1, 3) - vgt(3, 1))*eigvec(2) &
+                                 + (vgt(2, 1) - vgt(1, 2))*eigvec(3)
+
+                    ! As eigenvector can have +/- signs, we can choose the sign
+                    ! so that omega_proj is positive
+                    if (omega_proj < 0._wp) then
+                        eigvec = -eigvec
+                        omega_proj = -omega_proj
+                    end if
+
+                    ! Find imaginary part of complex eigenvalue
+                    lci = li(mod(idx, 3) + 1)
+
+                    ! Compute Liutex magnitude
+                    alpha = omega_proj**2._wp - 4._wp*lci**2._wp ! (2*alpha)^2
+                    if (alpha > 0._wp) then
+                        liutex_mag(j, k, l) = omega_proj - sqrt(alpha)
+                    else
+                        liutex_mag(j, k, l) = omega_proj
+                    end if
+
+                    ! Compute Liutex axis
+                    liutex_axis(j, k, l, 1) = eigvec(1)
+                    liutex_axis(j, k, l, 2) = eigvec(2)
+                    liutex_axis(j, k, l, 3) = eigvec(3)
+
+                end do
+            end do
+        end do
+
+    end subroutine s_derive_liutex
+
     !>  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_derived_variables.fpp

@@ -252,6 +252,7 @@ module m_global_parameters
     logical :: c_wrt
     logical, dimension(3) :: omega_wrt
     logical :: qm_wrt
+    logical :: liutex_wrt
     logical :: schlieren_wrt
     logical :: cf_wrt
     logical :: ib
@@ -432,6 +433,7 @@ contains
         c_wrt = .false.
         omega_wrt = .false.
         qm_wrt = .false.
+        liutex_wrt = .false.
         schlieren_wrt = .false.
         sim_data = .false.
         cf_wrt = .false.
@@ -835,7 +837,7 @@ contains
         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 .or. qm_wrt) then
+        if (any(omega_wrt) .or. schlieren_wrt .or. qm_wrt .or. liutex_wrt) then
             fd_number = max(1, fd_order/2)
             buff_size = buff_size + fd_number
         end if

src/post_process/m_global_parameters.fpp

@@ -101,7 +101,7 @@ contains
             & 'heat_ratio_wrt', 'pi_inf_wrt', 'pres_inf_wrt', 'cons_vars_wrt', &
             & 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt','chem_wrt_T', &
             & 'bubbles_euler', 'qbmm', 'polytropic', 'polydisperse',           &
-            & 'file_per_process', 'relax', 'cf_wrt', 'igr',                    &
+            & 'file_per_process', 'relax', 'cf_wrt', 'igr', 'liutex_wrt',      &
             & 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt',          &
             & 'surface_tension', 'hyperelasticity', 'bubbles_lagrange',        &
             & 'output_partial_domain', 'relativity', 'cont_damage', 'bc_io',   &