entropy_viscosity.f90

Go to the documentation of this file.

! Copyright (c) 2025-2026, 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 entropy_viscosity
  use num_types, only : rp
  use case, only : case_t
  use utils, only : neko_error
  use registry, only : neko_registry
  use artificial_viscosity_model, only : avm_t
  use json_module, only : json_file
  use json_utils, only : json_get_or_default
  use field, only : field_t
  use field_math, only : field_cfill, field_glsum, field_cadd, field_copy
  use field_series, only : field_series_t
  use coefs, only : coef_t
  use dofmap, only : dofmap_t
  use gather_scatter, only : gs_t
  use time_state, only : time_state_t
  use bdf_time_scheme, only : bdf_time_scheme_t
  use operators, only : div
  use math, only : glmax, absval, col2
  use scratch_registry, only : neko_scratch_registry
  use mesh, only : mesh_t
  use space, only : space_t
  use gather_scatter, only : gs_t
  use gs_ops, only : gs_op_add
  use neko_config, only : neko_bcknd_device
  use device, only : device_memcpy, host_to_device, device_to_host
  use device_math, only : device_col3, device_absval, device_glsum, device_col2
  use entropy_viscosity_cpu, only : entropy_viscosity_compute_residual_cpu, &
       entropy_viscosity_compute_viscosity_cpu, &
       entropy_viscosity_apply_element_max_cpu, &
       entropy_viscosity_clamp_to_low_order_cpu, &
       entropy_viscosity_smooth_divide_cpu
  use entropy_viscosity_device, only : &
       entropy_viscosity_compute_residual_device, &
       entropy_viscosity_compute_viscosity_device, &
       entropy_viscosity_apply_element_max_device, &
       entropy_viscosity_clamp_to_low_order_device, &
       entropy_viscosity_smooth_divide_device
  use compressible_ops_cpu, only : compressible_ops_cpu_compute_entropy
  use compressible_ops_device, only : compressible_ops_device_compute_entropy
  use fluid_scheme_compressible, only : fluid_scheme_compressible_t
  implicit none
  private
 
  type, public, extends(avm_t) :: entropy_viscosity_t
     real(kind=rp) :: c_avisc_entropy
     real(kind=rp) :: c_avisc_low
     real(kind=rp) :: gamma
     type(field_t) :: entropy_residual
     type(field_series_t) :: s_lag
     type(field_t), pointer :: s => null()
     type(field_t), pointer :: p => null()
     type(field_t), pointer :: rho => null()
     type(field_t), pointer :: u => null()
     type(field_t), pointer :: v => null()
     type(field_t), pointer :: w => null()
     type(field_t) :: h
     type(field_t), pointer :: max_wave_speed => null()
     type(mesh_t), pointer :: msh => null()
     type(space_t), pointer :: xh => null()
     type(gs_t), pointer :: gs => null()
   contains
     procedure, pass(this) :: init => entropy_viscosity_init
     procedure, pass(this) :: free => entropy_viscosity_free
     procedure, pass(this) :: preprocess => entropy_viscosity_preprocess
     procedure, pass(this) :: compute => entropy_viscosity_update_lag
     procedure, pass(this) :: restart => entropy_viscosity_restart
     procedure, pass(this) :: compute_h => entropy_viscosity_compute_h
     procedure, pass(this) :: set_fields => entropy_viscosity_set_fields
     procedure, pass(this) :: compute_entropy => &
          entropy_viscosity_compute_entropy
     procedure, pass(this), private :: compute_residual => &
          entropy_viscosity_compute_residual
     procedure, pass(this), private :: compute_viscosity => &
          entropy_viscosity_compute_viscosity
     procedure, pass(this), private :: smooth_viscosity => &
          entropy_viscosity_smooth_viscosity
     procedure, pass(this), private :: apply_element_max => &
          entropy_viscosity_apply_element_max
     procedure, pass(this), private :: low_order_viscosity => &
          entropy_viscosity_low_order
  end type entropy_viscosity_t
 
contains
 
  subroutine entropy_viscosity_init(this, case, json)
    class(entropy_viscosity_t), intent(inout) :: this
    type(json_file), intent(inout) :: json
    class(case_t), intent(inout), target :: case
    character(len=:), allocatable :: reg_coeff_name
 
    call this%free()
 
    select type (fluid => case%fluid)
    class is (fluid_scheme_compressible_t)
       call this%set_fields(fluid%p, fluid%rho, fluid%u, fluid%v, &
            fluid%w, fluid%max_wave_speed, fluid%msh, fluid%Xh, fluid%gs_Xh, &
            fluid%gamma)
    class default
       call neko_error('Entropy viscosity requires a compressible fluid scheme')
    end select
 
    call json_get_or_default(json, 'field_name', &
         reg_coeff_name, "entropy_viscosity")
    call this%init_base(case%fluid%dm_Xh, case%fluid%c_Xh, trim(reg_coeff_name))
    if (allocated(reg_coeff_name)) deallocate(reg_coeff_name)
 
    call json_get_or_default(json, 'c_avisc_low', this%c_avisc_low, 0.5_rp)
    call json_get_or_default(json, 'c_avisc_entropy', &
         this%c_avisc_entropy, 1.0_rp)
 
    call this%entropy_residual%init(this%dof, 'entropy_residual')
 
    call neko_registry%add_field(this%dof, 'S', .true.)
    this%S => neko_registry%get_field('S')
    call this%compute_entropy()
    call this%S_lag%init(this%S, 3)
 
    call this%h%init(this%dof, 'h')
    call this%compute_h()
 
  end subroutine entropy_viscosity_init
 
  subroutine entropy_viscosity_free(this)
    class(entropy_viscosity_t), intent(inout) :: this
 
    call this%free_base()
    call this%entropy_residual%free()
    call this%S_lag%free()
    call this%h%free()
 
    nullify(this%S)
    nullify(this%p)
    nullify(this%rho)
    nullify(this%u)
    nullify(this%v)
    nullify(this%w)
    nullify(this%max_wave_speed)
    nullify(this%msh)
    nullify(this%Xh)
    nullify(this%gs)
 
  end subroutine entropy_viscosity_free
 
  subroutine entropy_viscosity_preprocess(this, time)
    class(entropy_viscosity_t), intent(inout) :: this
    type(time_state_t), intent(in) :: time
 
    call this%compute_residual(time%tstep, &
         real(time%dt, kind=rp), &
         real(time%dtlag, kind=rp))
 
    call this%compute_viscosity(time%tstep)
 
  end subroutine entropy_viscosity_preprocess
 
  subroutine entropy_viscosity_compute_residual(this, tstep, dt, dt_lag)
    class(entropy_viscosity_t), intent(inout) :: this
    integer, intent(in) :: tstep
    real(kind=rp), intent(in) :: dt
    real(kind=rp), intent(in) :: dt_lag(10)
    integer :: n
    type(field_t), pointer :: us_field, vs_field, ws_field, div_field
    integer :: temp_indices(4)
    real(kind=rp) :: bdf_coeffs(4)
    type(bdf_time_scheme_t) :: bdf_scheme
    real(kind=rp) :: dt_local(10)
 
    if (tstep .le. 3) then
       return
    end if
 
    n = this%dof%size()
    call field_cfill(this%entropy_residual, 0.0_rp, n)
 
    bdf_coeffs = 0.0_rp
    dt_local = dt_lag
 
    call bdf_scheme%compute_coeffs(bdf_coeffs, dt_local, 3)
 
    if (neko_bcknd_device .eq. 1) then
       call entropy_viscosity_compute_residual_device( &
            this%entropy_residual%x_d, &
            this%S%x_d, this%S_lag%lf(1)%x_d, &
            this%S_lag%lf(2)%x_d, this%S_lag%lf(3)%x_d, &
            bdf_coeffs, dt, n)
    else
       call entropy_viscosity_compute_residual_cpu( &
            this%entropy_residual%x, &
            this%S%x, this%S_lag%lf(1)%x, &
            this%S_lag%lf(2)%x, this%S_lag%lf(3)%x, &
            bdf_coeffs, dt, n)
    end if
 
    call neko_scratch_registry%request_field(us_field, temp_indices(1), .false.)
    call neko_scratch_registry%request_field(vs_field, temp_indices(2), .false.)
    call neko_scratch_registry%request_field(ws_field, temp_indices(3), .false.)
    call neko_scratch_registry%request_field(div_field, temp_indices(4), &
         .false.)
 
    if (neko_bcknd_device .eq. 1) then
       call device_col3(us_field%x_d, this%u%x_d, this%S%x_d, n)
       call device_col3(vs_field%x_d, this%v%x_d, this%S%x_d, n)
       call device_col3(ws_field%x_d, this%w%x_d, this%S%x_d, n)
    else
       call entropy_viscosity_col3_vector_cpu(us_field%x, vs_field%x, &
            ws_field%x, this%u%x, this%v%x, this%w%x, this%S%x, n)
    end if
 
    call div(div_field%x, us_field%x, vs_field%x, ws_field%x, this%coef)
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%entropy_residual%x, this%entropy_residual%x_d, &
            n, device_to_host, sync = .false.)
       call device_memcpy(div_field%x, div_field%x_d, n, device_to_host, &
            sync = .true.)
    end if
 
    call entropy_viscosity_abs_add_cpu(this%entropy_residual%x, div_field%x, n)
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%entropy_residual%x, this%entropy_residual%x_d, &
            n, host_to_device, sync = .false.)
    end if
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine entropy_viscosity_compute_residual
 
  subroutine entropy_viscosity_compute_viscosity(this, tstep)
    class(entropy_viscosity_t), intent(inout) :: this
    integer, intent(in) :: tstep
    integer :: n, temp_indices(1)
    real(kind=rp) :: s_mean, n_s
    type(field_t), pointer :: temp_field
 
    n = this%dof%size()
 
    if (tstep .le. 3) then
       call field_cfill(this%reg_coeff, 0.0_rp, n)
       return
    end if
 
    call neko_scratch_registry%request_field(temp_field, temp_indices(1), &
         .false.)
 
    call field_cfill(temp_field, 1.0_rp, n)
    s_mean = field_glsum(this%S, n) / field_glsum(temp_field, n)
 
    call field_copy(temp_field, this%S, n)
    call field_cadd(temp_field, -s_mean, n)
 
    if (neko_bcknd_device .eq. 1) then
       call device_absval(temp_field%x_d, n)
       call device_memcpy(temp_field%x, temp_field%x_d, n, device_to_host, &
            sync = .true.)
    else
       call absval(temp_field%x, n)
    end if
 
    ! Normalization factor n_S = |S-S_mean|_inf
    n_s = glmax(temp_field%x, n)
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
    if (n_s < 1.0e-12_rp) then
       n_s = 1.0e-12_rp
    end if
 
    ! entropy viscosity = c_avisc_entropy * h^2 * entropy_residual / n_S
    if (neko_bcknd_device .eq. 1) then
       call entropy_viscosity_compute_viscosity_device( &
            this%reg_coeff%x_d, this%entropy_residual%x_d, &
            this%h%x_d, this%c_avisc_entropy, n_s, n)
    else
       call entropy_viscosity_compute_viscosity_cpu( &
            this%reg_coeff%x, this%entropy_residual%x, &
            this%h%x, this%c_avisc_entropy, n_s, n)
    end if
 
    ! artificial viscosity = min(entropy viscosity, low-order viscosity)
    if (neko_bcknd_device .eq. 1) then
       call entropy_viscosity_clamp_to_low_order_device( &
            this%reg_coeff%x_d, this%h%x_d, this%max_wave_speed%x_d, &
            this%c_avisc_low, n)
    else
       call entropy_viscosity_clamp_to_low_order_cpu( &
            this%reg_coeff%x, this%h%x, this%max_wave_speed%x, &
            this%c_avisc_low, n)
    end if
 
    call this%apply_element_max()
 
    call this%smooth_viscosity()
 
  end subroutine entropy_viscosity_compute_viscosity
 
  subroutine entropy_viscosity_smooth_viscosity(this)
    class(entropy_viscosity_t), intent(inout) :: this
    integer :: n
    type(field_t), pointer :: temp_field, mult_field
    integer :: temp_indices(2)
 
    n = this%dof%size()
 
    call neko_scratch_registry%request_field(temp_field, temp_indices(1), &
         .false.)
    call neko_scratch_registry%request_field(mult_field, temp_indices(2), &
         .false.)
 
    call field_copy(temp_field, this%reg_coeff, n)
    call this%gs%op(temp_field, gs_op_add)
 
    call field_cfill(mult_field, 1.0_rp, n)
    call this%gs%op(mult_field, gs_op_add)
 
    if (neko_bcknd_device .eq. 1) then
       call entropy_viscosity_smooth_divide_device( &
            this%reg_coeff%x_d, temp_field%x_d, mult_field%x_d, n)
    else
       call entropy_viscosity_smooth_divide_cpu( &
            this%reg_coeff%x, temp_field%x, mult_field%x, n)
    end if
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine entropy_viscosity_smooth_viscosity
 
  subroutine entropy_viscosity_apply_element_max(this)
    class(entropy_viscosity_t), intent(inout) :: this
    integer :: lx
 
    lx = this%Xh%lx
 
    if (neko_bcknd_device .eq. 1) then
       call entropy_viscosity_apply_element_max_device( &
            this%reg_coeff%x_d, lx, this%msh%nelv)
    else
       call entropy_viscosity_apply_element_max_cpu( &
            this%reg_coeff%x, lx, this%msh%nelv)
    end if
 
  end subroutine entropy_viscosity_apply_element_max
 
  subroutine entropy_viscosity_set_fields(this, p, rho, u, v, w, &
       max_wave_speed, msh, Xh, gs, gamma)
    class(entropy_viscosity_t), intent(inout) :: this
    type(field_t), target, intent(in) :: p, rho, u, v, w, max_wave_speed
    type(mesh_t), target, intent(in) :: msh
    type(space_t), target, intent(in) :: Xh
    type(gs_t), target, intent(in) :: gs
    real(kind=rp), intent(in) :: gamma
 
    this%p => p
    this%rho => rho
    this%u => u
    this%v => v
    this%w => w
    this%max_wave_speed => max_wave_speed
    this%msh => msh
    this%Xh => xh
    this%gs => gs
    this%gamma = gamma
 
  end subroutine entropy_viscosity_set_fields
 
  subroutine entropy_viscosity_compute_entropy(this)
    class(entropy_viscosity_t), intent(inout) :: this
    integer :: n
 
    n = this%dof%size()
 
    if (neko_bcknd_device .eq. 1) then
       call compressible_ops_device_compute_entropy(this%S, this%p, this%rho, &
            this%gamma, n)
    else
       call compressible_ops_cpu_compute_entropy(this%S%x, this%p%x, &
            this%rho%x, this%gamma, n)
    end if
 
  end subroutine entropy_viscosity_compute_entropy
 
  subroutine entropy_viscosity_restart(this, time)
    class(entropy_viscosity_t), intent(inout) :: this
    type(time_state_t), intent(in) :: time
 
    call this%compute_entropy()
    call this%S_lag%set(this%S)
 
  end subroutine entropy_viscosity_restart
 
  subroutine entropy_viscosity_update_lag(this, time)
    class(entropy_viscosity_t), intent(inout) :: this
    type(time_state_t), intent(in) :: time
 
    call this%S_lag%update()
    call this%compute_entropy()
 
 
  end subroutine entropy_viscosity_update_lag
 
  subroutine entropy_viscosity_col3_vector_cpu(us, vs, ws, u, v, w, S, n)
    integer, intent(in) :: n
    real(kind=rp), intent(out) :: us(n), vs(n), ws(n)
    real(kind=rp), intent(in) :: u(n), v(n), w(n), s(n)
    integer :: i
 
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do simd
    do i = 1, n
       us(i) = u(i) * s(i)
       vs(i) = v(i) * s(i)
       ws(i) = w(i) * s(i)
    end do
    !$omp end parallel do simd
  end subroutine entropy_viscosity_col3_vector_cpu
 
  subroutine entropy_viscosity_abs_add_cpu(entropy_residual, div_field, n)
    integer, intent(in) :: n
    real(kind=rp), intent(inout) :: entropy_residual(n)
    real(kind=rp), intent(in) :: div_field(n)
    integer :: i
 
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do simd
    do i = 1, n
       entropy_residual(i) = abs(entropy_residual(i) + div_field(i))
    end do
    !$omp end parallel do simd
  end subroutine entropy_viscosity_abs_add_cpu
 
  pure function entropy_viscosity_low_order(this, i) result(visc)
    class(entropy_viscosity_t), intent(in) :: this
    integer, intent(in) :: i
    real(kind=rp) :: visc
 
    visc = this%c_avisc_low * this%h%x(i,1,1,1) * this%max_wave_speed%x(i,1,1,1)
 
  end function entropy_viscosity_low_order
 
 
  subroutine entropy_viscosity_compute_h(this)
    class(entropy_viscosity_t), intent(inout) :: this
    integer :: e, i, j, k
    integer :: im, ip, jm, jp, km, kp
    real(kind=rp) :: di, dj, dk, ndim_inv
    integer :: lx_half, ly_half, lz_half
 
    lx_half = this%coef%Xh%lx / 2
    ly_half = this%coef%Xh%ly / 2
    lz_half = this%coef%Xh%lz / 2
 
    do concurrent(e = 1:this%coef%msh%nelv)
       do concurrent(k = 1:this%coef%Xh%lz, &
            j = 1:this%coef%Xh%ly, i = 1:this%coef%Xh%lx)
          km = max(1, k-1)
          kp = min(this%coef%Xh%lz, k+1)
 
          jm = max(1, j-1)
          jp = min(this%coef%Xh%ly, j+1)
 
          im = max(1, i-1)
          ip = min(this%coef%Xh%lx, i+1)
 
          di = (this%coef%dof%x%x(ip, j, k, e) - &
               this%coef%dof%x%x(im, j, k, e))**2 &
               + (this%coef%dof%y%x(ip, j, k, e) - &
               this%coef%dof%y%x(im, j, k, e))**2 &
               + (this%coef%dof%z%x(ip, j, k, e) - &
               this%coef%dof%z%x(im, j, k, e))**2
 
          dj = (this%coef%dof%x%x(i, jp, k, e) - &
               this%coef%dof%x%x(i, jm, k, e))**2 &
               + (this%coef%dof%y%x(i, jp, k, e) - &
               this%coef%dof%y%x(i, jm, k, e))**2 &
               + (this%coef%dof%z%x(i, jp, k, e) - &
               this%coef%dof%z%x(i, jm, k, e))**2
 
          dk = (this%coef%dof%x%x(i, j, kp, e) - &
               this%coef%dof%x%x(i, j, km, e))**2 &
               + (this%coef%dof%y%x(i, j, kp, e) - &
               this%coef%dof%y%x(i, j, km, e))**2 &
               + (this%coef%dof%z%x(i, j, kp, e) - &
               this%coef%dof%z%x(i, j, km, e))**2
 
          di = sqrt(di) / (ip - im)
          dj = sqrt(dj) / (jp - jm)
          dk = sqrt(dk) / (kp - km)
          this%h%x(i,j,k,e) = (di * dj * dk)**(1.0_rp / 3.0_rp)
 
       end do
    end do
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%h%x, this%h%x_d, this%h%dof%size(),&
            host_to_device, sync = .false.)
       call this%gs%op(this%h, gs_op_add)
       call device_col2(this%h%x_d, this%coef%mult_d, this%h%dof%size())
    else
       call this%gs%op(this%h, gs_op_add)
       call col2(this%h%x, this%coef%mult, this%h%dof%size())
    end if
 
  end subroutine entropy_viscosity_compute_h
 
end module entropy_viscosity