fluid_scheme_compressible.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 fluid_scheme_compressible
  use field, only : field_t
  use field_math, only : field_col2, field_col3, &
       field_cmult2, field_cmult, field_addcol3, field_add2, &
       field_cfill
 
  use registry, only : neko_registry
  use bc, only : bc_t
  use fluid_scheme_base, only : fluid_scheme_base_t
  use json_module, only : json_file
  use num_types, only : rp, dp
  use mesh, only : mesh_t
  use scratch_registry, only : neko_scratch_registry
  use space, only : gll
  use user_intf, only : user_t, user_material_properties_intf, &
       dummy_user_material_properties
  use json_utils, only : json_get_or_default, json_get_or_lookup, &
       json_get_or_lookup_or_default
  use mpi_f08
  use operators, only : cfl_compressible
  use device, only : device_memcpy, host_to_device
  use compressible_ops_cpu, only : &
       compressible_ops_cpu_compute_max_wave_speed
  use compressible_ops_device, only : &
       compressible_ops_device_compute_max_wave_speed
  use neko_config, only : neko_bcknd_device
  use time_state, only : time_state_t
  use logger, only : neko_log, log_size
  use math, only : glsum, glamax
  use device_math, only : device_glamax
  use compressible_residual, only : compressible_rhs_set_physical_flux
  use num_types, only : i8
  implicit none
  private
 
  type, public, abstract, extends(fluid_scheme_base_t) :: &
       fluid_scheme_compressible_t
     type(field_t), pointer :: m_x => null() 
     type(field_t), pointer :: m_y => null() 
     type(field_t), pointer :: m_z => null() 
     type(field_t), pointer :: e => null() 
     type(field_t), pointer :: temperature => null() 
     type(field_t), pointer :: max_wave_speed => null() 
     type(field_t), pointer :: artificial_visc => null() 
     type(field_t), pointer :: kappa => null() 
 
     real(kind=rp) :: gamma
 
     real(kind=rp) :: mu_amax = 0.0_rp
     real(kind=rp) :: kappa_amax = 0.0_rp
 
     logical :: add_physical_flux = .false.
     logical :: add_physical_stress = .false.
     logical :: variable_material_properties = .false.
 
     integer(kind=i8) :: glb_n_points
     integer(kind=i8) :: glb_unique_points
 
   contains
     procedure, pass(this) :: scheme_init => fluid_scheme_compressible_init
     procedure, pass(this) :: scheme_free => fluid_scheme_compressible_free
 
     procedure, pass(this) :: validate => fluid_scheme_compressible_validate
     procedure, pass(this) :: compute_cfl &
          => fluid_scheme_compressible_compute_cfl
     procedure, pass(this) :: set_material_properties => &
          fluid_scheme_compressible_set_material_properties
     procedure, pass(this) :: update_material_properties => &
          fluid_scheme_compressible_update_material_properties
     procedure, pass(this) :: update_physical_flux => &
          fluid_scheme_compressible_update_physical_flux
     procedure, pass(this) :: compute_max_wave_speed => &
          fluid_scheme_compressible_compute_max_wave_speed
     procedure, pass(this) :: log_solver_info => &
          fluid_scheme_compressible_log_solver_info
 
  end type fluid_scheme_compressible_t
 
contains
  subroutine fluid_scheme_compressible_init(this, msh, lx, params, scheme, user)
    class(fluid_scheme_compressible_t), target, intent(inout) :: this
    type(mesh_t), target, intent(inout) :: msh
    integer, intent(in) :: lx
    character(len=*), intent(in) :: scheme
    type(json_file), target, intent(inout) :: params
    type(user_t), target, intent(in) :: user
    logical :: compress_geo
    real(kind=rp) :: compress_tol
 
    !
    ! SEM simulation fundamentals
    !
 
    this%msh => msh
 
    if (msh%gdim .eq. 2) then
       call this%Xh%init(gll, lx, lx)
    else
       call this%Xh%init(gll, lx, lx, lx)
    end if
 
    call this%dm_Xh%init(msh, this%Xh)
 
    call this%gs_Xh%init(this%dm_Xh)
 
    call this%c_Xh%init(this%gs_Xh)
 
    ! Assign Dofmap to scratch registry
    call neko_scratch_registry%set_dofmap(this%dm_Xh)
 
    ! Case parameters
    this%params => params
 
    ! Assign a name
    call json_get_or_default(params, 'case.fluid.name', this%name, "fluid")
 
    ! Material properties will be set up via set_material_properties
    call neko_registry%add_field(this%dm_Xh, this%name // "_rho")
    this%rho => neko_registry%get_field(this%name // "_rho")
 
    ! Assign momentum fields
    call neko_registry%add_field(this%dm_Xh, "m_x")
    call neko_registry%add_field(this%dm_Xh, "m_y")
    call neko_registry%add_field(this%dm_Xh, "m_z")
    this%m_x => neko_registry%get_field("m_x")
    this%m_y => neko_registry%get_field("m_y")
    this%m_z => neko_registry%get_field("m_z")
    call this%m_x%init(this%dm_Xh, "m_x")
    call this%m_y%init(this%dm_Xh, "m_y")
    call this%m_z%init(this%dm_Xh, "m_z")
 
    ! Assign energy field
    call neko_registry%add_field(this%dm_Xh, "E")
    this%E => neko_registry%get_field("E")
    call this%E%init(this%dm_Xh, "E")
 
    ! Assign temperature field
    call neko_registry%add_field(this%dm_Xh, "temperature")
    this%temperature => neko_registry%get_field("temperature")
    call this%temperature%init(this%dm_Xh, "temperature")
 
    ! Assign maximum wave speed field
    call neko_registry%add_field(this%dm_Xh, "max_wave_speed")
    this%max_wave_speed => neko_registry%get_field("max_wave_speed")
    call this%max_wave_speed%init(this%dm_Xh, "max_wave_speed")
 
    ! Assign artificial viscosity field (without physical viscosity)
    call neko_registry%add_field(this%dm_Xh, "artificial_visc")
    this%artificial_visc => neko_registry%get_field("artificial_visc")
    call this%artificial_visc%init(this%dm_Xh, "artificial_visc")
 
    ! ! Assign velocity fields
    call neko_registry%add_field(this%dm_Xh, "u")
    call neko_registry%add_field(this%dm_Xh, "v")
    call neko_registry%add_field(this%dm_Xh, "w")
    this%u => neko_registry%get_field("u")
    this%v => neko_registry%get_field("v")
    this%w => neko_registry%get_field("w")
    call this%u%init(this%dm_Xh, "u")
    call this%v%init(this%dm_Xh, "v")
    call this%w%init(this%dm_Xh, "w")
    call neko_registry%add_field(this%dm_Xh, 'p')
    this%p => neko_registry%get_field('p')
    call this%p%init(this%dm_Xh, "p")
 
    !
    ! Setup right-hand side fields.
    !
    allocate(this%f_x)
    allocate(this%f_y)
    allocate(this%f_z)
    call this%f_x%init(this%dm_Xh, fld_name = "fluid_rhs_x")
    call this%f_y%init(this%dm_Xh, fld_name = "fluid_rhs_y")
    call this%f_z%init(this%dm_Xh, fld_name = "fluid_rhs_z")
 
    ! Material properties
    call this%set_material_properties(params, user)
 
    ! Compressible parameters
    call json_get_or_default(params, 'case.fluid.gamma', this%gamma, 1.4_rp)
 
    ! Calculate global points for logging
    this%glb_n_points = int(this%msh%glb_nelv, i8)*int(this%Xh%lxyz, i8)
    this%glb_unique_points = int(glsum(this%c_Xh%mult, this%dm_Xh%size()), i8)
 
    !
    ! Log solver information
    !
    call this%log_solver_info(params, scheme, lx)
 
    ! Compressed geometric factors, read by the standard Helmholtz operator
    call json_get_or_default(params, &
         'case.numerics.compress_geometric_factors', compress_geo, .false.)
    if (compress_geo) then
       if (params%valid_path( &
            'case.numerics.compress_geometric_factors_tolerance')) then
          call json_get_or_lookup(params, &
               'case.numerics.compress_geometric_factors_tolerance', &
               compress_tol)
          call this%c_Xh%enable_geo_compression(compress_tol)
       else
          call this%c_Xh%enable_geo_compression()
       end if
    end if
  end subroutine fluid_scheme_compressible_init
 
  subroutine fluid_scheme_compressible_free(this)
    class(fluid_scheme_compressible_t), intent(inout) :: this
    class(bc_t), pointer :: bc
    integer :: i
 
    do i = 1, this%bcs_vel%size()
       bc => this%bcs_vel%get(i)
       if (associated(bc)) then
          call bc%free()
          deallocate(bc)
       end if
    end do
    call this%bcs_vel%free()
 
    do i = 1, this%bcs_prs%size()
       bc => this%bcs_prs%get(i)
       if (associated(bc)) then
          call bc%free()
          deallocate(bc)
       end if
    end do
    call this%bcs_prs%free()
 
    call this%dm_Xh%free()
    call this%gs_Xh%free()
    call this%c_Xh%free()
    call this%Xh%free()
 
    if (associated(this%m_x)) then
       call this%m_x%free()
    end if
 
    if (associated(this%m_y)) then
       call this%m_y%free()
    end if
 
    if (associated(this%m_z)) then
       call this%m_z%free()
    end if
 
    if (associated(this%E)) then
       call this%E%free()
    end if
 
    if (associated(this%temperature)) then
       call this%temperature%free()
    end if
 
    if (associated(this%max_wave_speed)) then
       call this%max_wave_speed%free()
    end if
 
    if (associated(this%f_x)) then
       call this%f_x%free()
       deallocate(this%f_x)
    end if
 
    if (associated(this%f_y)) then
       call this%f_y%free()
       deallocate(this%f_y)
    end if
 
    if (associated(this%f_z)) then
       call this%f_z%free()
       deallocate(this%f_z)
    end if
 
    nullify(this%f_x)
    nullify(this%f_y)
    nullify(this%f_z)
 
    nullify(this%m_x)
    nullify(this%m_y)
    nullify(this%m_z)
    nullify(this%E)
    nullify(this%temperature)
    nullify(this%max_wave_speed)
 
    nullify(this%u)
    nullify(this%v)
    nullify(this%w)
    nullify(this%p)
    nullify(this%rho)
    nullify(this%mu)
    nullify(this%kappa)
 
    call this%material_properties%free()
 
  end subroutine fluid_scheme_compressible_free
 
  subroutine fluid_scheme_compressible_validate(this)
    class(fluid_scheme_compressible_t), target, intent(inout) :: this
    integer :: n, i
    type(field_t), pointer :: temp
    integer :: temp_indices(1)
 
    n = this%dm_Xh%size()
    call neko_scratch_registry%request_field(temp, temp_indices(1), .false.)
 
    call field_col3(this%m_x, this%u, this%rho)
    call field_col3(this%m_y, this%v, this%rho)
    call field_col3(this%m_z, this%w, this%rho)
 
    call field_cmult2(this%E, this%p, 1.0_rp/(this%gamma - 1.0_rp), n)
    call field_col3(temp, this%u, this%u, n)
    call field_addcol3(temp, this%v, this%v, n)
    call field_addcol3(temp, this%w, this%w, n)
    call field_col2(temp, this%rho, n)
    call field_cmult(temp, 0.5_rp, n)
    call field_add2(this%E, temp, n)
 
    do i = 1, n
       this%temperature%x(i,1,1,1) = this%p%x(i,1,1,1) / &
            (this%rho%x(i,1,1,1) * (this%gamma - 1.0_rp))
    end do
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
    call this%compute_max_wave_speed()
 
  end subroutine fluid_scheme_compressible_validate
 
  function fluid_scheme_compressible_compute_cfl(this, dt) result(c)
    class(fluid_scheme_compressible_t), intent(in) :: this
    real(kind=dp), intent(in) :: dt
    real(kind=dp) :: c
    integer :: n
 
    associate(u => this%u, v => this%v, w => this%w, p => this%p, &
         rho => this%rho, xh => this%Xh, c_xh => this%c_Xh, &
         msh => this%msh, gamma => this%gamma, &
         max_wave_speed => this%max_wave_speed)
 
      n = xh%lx * xh%ly * xh%lz * msh%nelv
 
      ! Use the compressible CFL function with precomputed maximum wave speed
      c = cfl_compressible(dt, max_wave_speed%x, xh, c_xh, msh%nelv, msh%gdim)
    end associate
 
  end function fluid_scheme_compressible_compute_cfl
 
  subroutine fluid_scheme_compressible_set_material_properties(this, &
       params, user)
    class(fluid_scheme_compressible_t), target, intent(inout) :: this
    type(json_file), intent(inout) :: params
    type(user_t), target, intent(in) :: user
    procedure(user_material_properties_intf), pointer :: dummy_mp_ptr
    type(time_state_t) :: dummy_time_state
    real(kind=rp) :: const_mu, const_kappa
 
    dummy_mp_ptr => dummy_user_material_properties
 
    call neko_registry%add_field(this%dm_Xh, this%name // "_mu")
    this%mu => neko_registry%get_field(this%name // "_mu")
 
    call neko_registry%add_field(this%dm_Xh, this%name // "_kappa")
    this%kappa => neko_registry%get_field(this%name // "_kappa")
 
    call this%material_properties%init(3)
    call this%material_properties%assign(1, this%rho)
    call this%material_properties%assign(2, this%mu)
    call this%material_properties%assign(3, this%kappa)
 
    if (.not. associated(user%material_properties, dummy_mp_ptr)) then
       this%user_material_properties => user%material_properties
       this%variable_material_properties = .true.
       call user%material_properties(this%name, this%material_properties, &
            dummy_time_state)
    else
       this%variable_material_properties = .false.
       this%user_material_properties => dummy_user_material_properties
       call json_get_or_lookup_or_default(params, 'case.fluid.mu', const_mu, &
            0.0_rp)
       call json_get_or_lookup_or_default(params, 'case.fluid.kappa', &
            const_kappa, 0.0_rp)
 
       call field_cfill(this%mu, const_mu)
       call field_cfill(this%kappa, const_kappa)
    end if
 
    ! Reported by log_solver_info, once the Fluid section is open.
    call this%update_physical_flux()
 
  end subroutine fluid_scheme_compressible_set_material_properties
 
  subroutine fluid_scheme_compressible_update_material_properties(this, time)
    use device, only : device_memcpy, host_to_device
    class(fluid_scheme_compressible_t), intent(inout) :: this
    type(time_state_t), intent(in) :: time
 
    if (.not. this%variable_material_properties) return
 
    call this%user_material_properties(this%name, this%material_properties, &
         time)
 
    call this%update_physical_flux()
 
  end subroutine fluid_scheme_compressible_update_material_properties
 
  subroutine fluid_scheme_compressible_update_physical_flux(this)
    class(fluid_scheme_compressible_t), intent(inout) :: this
    real(kind=rp) :: mu_amax, kappa_amax
    integer :: n
 
    n = this%dm_Xh%size()
 
    if (neko_bcknd_device .eq. 1) then
       mu_amax = device_glamax(this%mu%x_d, n)
       kappa_amax = device_glamax(this%kappa%x_d, n)
    else
       mu_amax = glamax(this%mu%x, n)
       kappa_amax = glamax(this%kappa%x, n)
    end if
 
    this%mu_amax = mu_amax
    this%kappa_amax = kappa_amax
 
    this%add_physical_stress = (mu_amax .ne. 0.0_rp)
    this%add_physical_flux = this%add_physical_stress .or. &
         (kappa_amax .ne. 0.0_rp)
 
    call compressible_rhs_set_physical_flux(this%add_physical_flux, &
         this%add_physical_stress)
 
  end subroutine fluid_scheme_compressible_update_physical_flux
 
  subroutine fluid_scheme_compressible_compute_max_wave_speed(this)
    class(fluid_scheme_compressible_t), intent(inout) :: this
    integer :: n
 
    n = this%u%dof%size()
 
    if (neko_bcknd_device .eq. 1) then
       call compressible_ops_device_compute_max_wave_speed( &
            this%max_wave_speed, this%u, this%v, this%w, this%gamma, this%p, &
            this%rho, n)
    else
       call compressible_ops_cpu_compute_max_wave_speed(this%max_wave_speed%x, &
            this%u%x, this%v%x, this%w%x, this%gamma, this%p%x, this%rho%x, n)
    end if
 
  end subroutine fluid_scheme_compressible_compute_max_wave_speed
 
  subroutine fluid_scheme_compressible_log_solver_info(this, params, scheme, lx)
    class(fluid_scheme_compressible_t), intent(inout) :: this
    type(json_file), intent(inout) :: params
    character(len=*), intent(in) :: scheme
    integer, intent(in) :: lx
    character(len=LOG_SIZE) :: log_buf
    logical :: logical_val
    integer :: integer_val
 
    call neko_log%section('Fluid')
    write(log_buf, '(A, A)') 'Type       : ', trim(scheme)
    call neko_log%message(log_buf)
    write(log_buf, '(A, A)') 'Name       : ', trim(this%name)
    call neko_log%message(log_buf)
 
    ! Polynomial order
    if (lx .lt. 10) then
       write(log_buf, '(A, I1)') 'Poly order : ', lx-1
    else if (lx .ge. 10) then
       write(log_buf, '(A, I2)') 'Poly order : ', lx-1
    else
       write(log_buf, '(A, I3)') 'Poly order : ', lx-1
    end if
    call neko_log%message(log_buf)
 
    ! Global points information
    write(log_buf, '(A, I0)') 'GLL points : ', this%glb_n_points
    call neko_log%message(log_buf)
    write(log_buf, '(A, I0)') 'Unique pts.: ', this%glb_unique_points
    call neko_log%message(log_buf)
 
    ! Material properties
    write(log_buf, '(A,ES13.6)') 'gamma      :', this%gamma
    call neko_log%message(log_buf)
 
    call json_get_or_default(params, 'case.numerics.time_order', integer_val, 4)
    write(log_buf, '(A, I0)') 'RK order   : ', integer_val
    call neko_log%message(log_buf)
 
    call json_get_or_default(params, &
         'case.numerics.compress_geometric_factors', logical_val, .false.)
    write(log_buf, '(A, L1)') 'Compress G : ', logical_val
    call neko_log%message(log_buf)
 
    ! Physical viscosity and conductivity. For user-defined properties the
    ! fields need not be uniform, so report the largest magnitude found.
    if (this%variable_material_properties) then
       write(log_buf, '(A,ES13.6,A)') 'mu         :', this%mu_amax, &
            ' (user, max)'
       call neko_log%message(log_buf)
       write(log_buf, '(A,ES13.6,A)') 'kappa      :', this%kappa_amax, &
            ' (user, max)'
       call neko_log%message(log_buf)
    else
       write(log_buf, '(A,ES13.6)') 'mu         :', this%mu_amax
       call neko_log%message(log_buf)
       write(log_buf, '(A,ES13.6)') 'kappa      :', this%kappa_amax
       call neko_log%message(log_buf)
    end if
 
    if (this%add_physical_flux) then
       write(log_buf, '(A, A)') 'NS fluxes  : ', 'enabled'
    else
       write(log_buf, '(A, A)') 'NS fluxes  : ', &
            'disabled (mu = kappa = 0, Euler)'
    end if
    call neko_log%message(log_buf)
 
    call neko_log%end_section()
 
  end subroutine fluid_scheme_compressible_log_solver_info
 
end module fluid_scheme_compressible