dynamic_smagorinsky_cpu.f90

Go to the documentation of this file.

! Copyright (c) 2024, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module dynamic_smagorinsky_cpu
  use num_types, only : rp
  use field_list, only : field_list_t
  use math, only : cadd, neko_eps, col2, sub2, col3, cmult
  use scratch_registry, only : neko_scratch_registry
  use registry, only : neko_registry
  use field, only : field_t
  use operators, only : strain_rate
  use coefs, only : coef_t
  use elementwise_filter, only : elementwise_filter_t
  use gs_ops, only : gs_op_add
  implicit none
  private
 
  public :: dynamic_smagorinsky_compute_cpu
 
contains
 
  subroutine dynamic_smagorinsky_compute_cpu(if_ext, t, tstep, coef, nut, delta, &
       c_dyn, test_filter, mij, lij, num, den)
    logical, intent(in) :: if_ext
    real(kind=rp), intent(in) :: t
    integer, intent(in) :: tstep
    type(coef_t), intent(in) :: coef
    type(field_t), intent(inout) :: nut
    type(field_t), intent(in) :: delta
    type(field_t), intent(inout) :: c_dyn
    type(elementwise_filter_t), intent(inout) :: test_filter
    type(field_t), intent(inout) :: mij(6), lij(6)
    type(field_t), intent(inout) :: num, den
 
    type(field_t), pointer :: u, v, w
    type(field_t) :: c_dyn_curr
    ! the strain rate tensor
    type(field_t), pointer :: s11, s22, s33, s12, s13, s23, s_abs
    real(kind=rp) :: alpha ! running averaging coefficient
    integer :: temp_indices(7)
    integer :: i
 
    if (tstep .eq. 1) then
       alpha = 1.0_rp
    else
       alpha = 0.9_rp
    end if
 
    if (if_ext .eqv. .true.) then
       u => neko_registry%get_field_by_name("u_e")
       v => neko_registry%get_field_by_name("v_e")
       w => neko_registry%get_field_by_name("w_e")
    else
       u => neko_registry%get_field_by_name("u")
       v => neko_registry%get_field_by_name("v")
       w => neko_registry%get_field_by_name("w")
    end if
 
 
    call neko_scratch_registry%request_field(s11, temp_indices(1), .false.)
    call neko_scratch_registry%request_field(s22, temp_indices(2), .false.)
    call neko_scratch_registry%request_field(s33, temp_indices(3), .false.)
    call neko_scratch_registry%request_field(s12, temp_indices(4), .false.)
    call neko_scratch_registry%request_field(s13, temp_indices(5), .false.)
    call neko_scratch_registry%request_field(s23, temp_indices(6), .false.)
    call neko_scratch_registry%request_field(s_abs, temp_indices(7), .false.)
 
    ! Compute the strain rate tensor
    call strain_rate(s11%x, s22%x, s33%x, s12%x, s13%x, s23%x, &
         u%x, v%x, w%x, coef)
 
    call coef%gs_h%op(s11%x, s11%dof%size(), gs_op_add)
    call coef%gs_h%op(s22%x, s11%dof%size(), gs_op_add)
    call coef%gs_h%op(s33%x, s11%dof%size(), gs_op_add)
    call coef%gs_h%op(s12%x, s11%dof%size(), gs_op_add)
    call coef%gs_h%op(s13%x, s11%dof%size(), gs_op_add)
    call coef%gs_h%op(s23%x, s11%dof%size(), gs_op_add)
 
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do
    do i = 1, u%dof%size()
       s11%x(i,1,1,1) = s11%x(i,1,1,1) * coef%mult(i,1,1,1)
       s22%x(i,1,1,1) = s22%x(i,1,1,1) * coef%mult(i,1,1,1)
       s33%x(i,1,1,1) = s33%x(i,1,1,1) * coef%mult(i,1,1,1)
       s12%x(i,1,1,1) = s12%x(i,1,1,1) * coef%mult(i,1,1,1)
       s13%x(i,1,1,1) = s13%x(i,1,1,1) * coef%mult(i,1,1,1)
       s23%x(i,1,1,1) = s23%x(i,1,1,1) * coef%mult(i,1,1,1)
    end do
    !$omp end parallel do
 
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do
    do i = 1, u%dof%size()
       s_abs%x(i,1,1,1) = sqrt(2.0_rp * (s11%x(i,1,1,1)*s11%x(i,1,1,1) + &
            s22%x(i,1,1,1)*s22%x(i,1,1,1) + &
            s33%x(i,1,1,1)*s33%x(i,1,1,1)) + &
            4.0_rp * (s12%x(i,1,1,1)*s12%x(i,1,1,1) + &
            s13%x(i,1,1,1)*s13%x(i,1,1,1) + &
            s23%x(i,1,1,1)*s23%x(i,1,1,1)))
    end do
    !$omp end parallel do
 
    call compute_lij_cpu(lij, u, v, w, test_filter, u%dof%size())
    call compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
         s_abs, test_filter, delta, u%dof%size())
    call compute_num_den_cpu(num, den, lij, mij, alpha, u%dof%size())
 
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do
    do i = 1, u%dof%size()
       if (den%x(i,1,1,1) .gt. 0.0_rp) then
          c_dyn%x(i,1,1,1) = 0.5_rp * (num%x(i,1,1,1)/den%x(i,1,1,1))
       else
          c_dyn%x(i,1,1,1) = 0.0_rp
       end if
       c_dyn%x(i,1,1,1) = max(c_dyn%x(i,1,1,1),0.0_rp)
       nut%x(i,1,1,1) = c_dyn%x(i,1,1,1) * delta%x(i,1,1,1)**2 &
            * s_abs%x(i,1,1,1)
    end do
    !$omp end parallel do
 
    call coef%gs_h%op(nut, gs_op_add)
    call col2(nut%x, coef%mult, nut%dof%size())
 
    call neko_scratch_registry%relinquish_field(temp_indices)
  end subroutine dynamic_smagorinsky_compute_cpu
 
  subroutine compute_lij_cpu(lij, u, v, w, test_filter, n)
    type(field_t), intent(inout) :: lij(6)
    type(field_t), pointer, intent(in) :: u, v, w
    type(elementwise_filter_t), intent(inout) :: test_filter
    integer, intent(in) :: n
    integer :: i
    integer :: temp_indices(3)
    type(field_t), pointer :: fu, fv, fw
 
    ! Use test filter for the velocity fields
    call neko_scratch_registry%request_field(fu, temp_indices(1), .false.)
    call neko_scratch_registry%request_field(fv, temp_indices(2), .false.)
    call neko_scratch_registry%request_field(fw, temp_indices(3), .false.)
    call test_filter%apply(fu, u)
    call test_filter%apply(fv, v)
    call test_filter%apply(fw, w)
 
    !! The first term
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do
    do i = 1, n
       lij(1)%x(i,1,1,1) = fu%x(i,1,1,1) * fu%x(i,1,1,1)
       lij(2)%x(i,1,1,1) = fv%x(i,1,1,1) * fv%x(i,1,1,1)
       lij(3)%x(i,1,1,1) = fw%x(i,1,1,1) * fw%x(i,1,1,1)
       lij(4)%x(i,1,1,1) = fu%x(i,1,1,1) * fv%x(i,1,1,1)
       lij(5)%x(i,1,1,1) = fu%x(i,1,1,1) * fw%x(i,1,1,1)
       lij(6)%x(i,1,1,1) = fv%x(i,1,1,1) * fw%x(i,1,1,1)
    end do
    !$omp end parallel do
 
    !! Subtract the second term:
    !! use test filter for the cross terms
    !! fu and fv are used as work array
    call col3(fu%x, u%x, u%x, n)
    call test_filter%apply(fv, fu)
    call sub2(lij(1)%x, fv%x, n)
 
    call col3(fu%x, v%x, v%x, n)
    call test_filter%apply(fv, fu)
    call sub2(lij(2)%x, fv%x, n)
 
    call col3(fu%x, w%x, w%x, n)
    call test_filter%apply(fv, fu)
    call sub2(lij(3)%x, fv%x, n)
 
    call col3(fu%x, u%x, v%x, n)
    call test_filter%apply(fv, fu)
    call sub2(lij(4)%x, fv%x, n)
 
    call col3(fu%x, u%x, w%x, n)
    call test_filter%apply(fv, fu)
    call sub2(lij(5)%x, fv%x, n)
 
    call col3(fu%x, v%x, w%x, n)
    call test_filter%apply(fv, fu)
    call sub2(lij(6)%x, fv%x, n)
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine compute_lij_cpu
 
  subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
       s_abs, test_filter, delta, n)
    type(field_t), intent(inout) :: mij(6)
    type(field_t), intent(inout) :: s11, s22, s33, s12, s13, s23, s_abs
    type(elementwise_filter_t), intent(inout) :: test_filter
    type(field_t), intent(in) :: delta
    integer, intent(in) :: n
 
    integer :: temp_indices(7)
    type(field_t), pointer :: fs11, fs22, fs33, fs12, fs13, fs23, fs_abs
    real(kind=rp) :: delta_ratio2 !test- to grid- filter ratio, squared
    integer :: i
    real(kind=rp) :: delta2
 
    delta_ratio2 = ((test_filter%nx-1.0_rp)/(test_filter%nt-1.0_rp))**2
 
    call neko_scratch_registry%request_field(fs11, temp_indices(1), .false.)
    call neko_scratch_registry%request_field(fs22, temp_indices(2), .false.)
    call neko_scratch_registry%request_field(fs33, temp_indices(3), .false.)
    call neko_scratch_registry%request_field(fs12, temp_indices(4), .false.)
    call neko_scratch_registry%request_field(fs13, temp_indices(5), .false.)
    call neko_scratch_registry%request_field(fs23, temp_indices(6), .false.)
    call neko_scratch_registry%request_field(fs_abs, temp_indices(7), .false.)
    !! The first term:
    !!                      _____ ____
    !! (delta_test/delta)^2 s_abs*s_ij
    call test_filter%apply(fs_abs, s_abs)
 
    call test_filter%apply(fs11, s11)
    call col3(mij(1)%x, fs_abs%x, fs11%x, n)
    call cmult(mij(1)%x, delta_ratio2, n)
 
    call test_filter%apply(fs22, s22)
    call col3(mij(2)%x, fs_abs%x, fs22%x, n)
    call cmult(mij(2)%x, delta_ratio2, n)
 
    call test_filter%apply(fs33, s33)
    call col3(mij(3)%x, fs_abs%x, fs33%x, n)
    call cmult(mij(3)%x, delta_ratio2, n)
 
    call test_filter%apply(fs12, s12)
    call col3(mij(4)%x, fs_abs%x, fs12%x, n)
    call cmult(mij(4)%x, delta_ratio2, n)
 
    call test_filter%apply(fs13, s13)
    call col3(mij(5)%x, fs_abs%x, fs13%x, n)
    call cmult(mij(5)%x, delta_ratio2, n)
 
    call test_filter%apply(fs23, s23)
    call col3(mij(6)%x, fs_abs%x, fs23%x, n)
    call cmult(mij(6)%x, delta_ratio2, n)
 
    !! Substract the second term:
    !!                      _____ ____   __________
    !! (delta_test/delta)^2 s_abs*s_ij - s_abs*s_ij
    !! fs11 and fs22 are used as work array
    call col3(fs11%x, s_abs%x, s11%x, n)
    call test_filter%apply(fs22, fs11)
    call sub2(mij(1)%x, fs22%x, n)
 
    call col3(fs11%x, s_abs%x, s22%x, n)
    call test_filter%apply(fs22, fs11)
    call sub2(mij(2)%x, fs22%x, n)
 
    call col3(fs11%x, s_abs%x, s33%x, n)
    call test_filter%apply(fs22, fs11)
    call sub2(mij(3)%x, fs22%x, n)
 
    call col3(fs11%x, s_abs%x, s12%x, n)
    call test_filter%apply(fs22, fs11)
    call sub2(mij(4)%x, fs22%x, n)
 
    call col3(fs11%x, s_abs%x, s13%x, n)
    call test_filter%apply(fs22, fs11)
    call sub2(mij(5)%x, fs22%x, n)
 
    call col3(fs11%x, s_abs%x, s23%x, n)
    call test_filter%apply(fs22, fs11)
    call sub2(mij(6)%x, fs22%x, n)
 
    !! Lastly multiplied by delta^2
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do private(i, delta2)
    do i = 1, n
       delta2 = delta%x(i,1,1,1)**2
       mij(1)%x(i,1,1,1) = mij(1)%x(i,1,1,1) * delta2
       mij(2)%x(i,1,1,1) = mij(2)%x(i,1,1,1) * delta2
       mij(3)%x(i,1,1,1) = mij(3)%x(i,1,1,1) * delta2
       mij(4)%x(i,1,1,1) = mij(4)%x(i,1,1,1) * delta2
       mij(5)%x(i,1,1,1) = mij(5)%x(i,1,1,1) * delta2
       mij(6)%x(i,1,1,1) = mij(6)%x(i,1,1,1) * delta2
    end do
    !$omp end parallel do
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine compute_mij_cpu
 
  subroutine compute_num_den_cpu(num, den, lij, mij, alpha, n)
    type(field_t), intent(inout) :: num, den
    type(field_t), intent(in) :: lij(6), mij(6)
    real(kind=rp), intent(in) :: alpha
    integer, intent(in) :: n
 
    real(kind=rp) :: num_curr, den_curr
    integer :: i
 
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do private(i, num_curr, den_curr)
    do i = 1, n
       num_curr = mij(1)%x(i,1,1,1)*lij(1)%x(i,1,1,1) + &
            mij(2)%x(i,1,1,1)*lij(2)%x(i,1,1,1) + &
            mij(3)%x(i,1,1,1)*lij(3)%x(i,1,1,1) + &
            2.0_rp*(mij(4)%x(i,1,1,1)*lij(4)%x(i,1,1,1) + &
            mij(5)%x(i,1,1,1)*lij(5)%x(i,1,1,1) + &
            mij(6)%x(i,1,1,1)*lij(6)%x(i,1,1,1))
       den_curr = mij(1)%x(i,1,1,1)*mij(1)%x(i,1,1,1) + &
            mij(2)%x(i,1,1,1)*mij(2)%x(i,1,1,1) + &
            mij(3)%x(i,1,1,1)*mij(3)%x(i,1,1,1) + &
            2.0_rp*(mij(4)%x(i,1,1,1)*mij(4)%x(i,1,1,1) + &
            mij(5)%x(i,1,1,1)*mij(5)%x(i,1,1,1) + &
            mij(6)%x(i,1,1,1)*mij(6)%x(i,1,1,1))
 
       ! running average over time
       num%x(i,1,1,1) = alpha * num%x(i,1,1,1) + (1.0_rp - alpha) * num_curr
       den%x(i,1,1,1) = alpha * den%x(i,1,1,1) + (1.0_rp - alpha) * den_curr
    end do
    !$omp end parallel do
 
  end subroutine compute_num_den_cpu
 
end module dynamic_smagorinsky_cpu