fluid_scheme_compressible_euler.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_euler
  use comm, only : neko_comm
  use advection, only : advection_t
  use device, only : device_memcpy, host_to_device
  use operators, only : div, rotate_cyc
  use field_series, only : field_series_t
  use bdf_time_scheme, only : bdf_time_scheme_t
  use time_scheme_controller, only : time_scheme_controller_t
  use math, only : col2
  use device_math, only : device_col2
  use field, only : field_t
  use fluid_scheme_compressible, only : fluid_scheme_compressible_t
  use scratch_registry, only : neko_scratch_registry
  use gs_ops, only : gs_op_add, gs_op_min, gs_op_max
  use gather_scatter, only : gs_t
  use num_types, only : rp
  use mesh, only : mesh_t
  use checkpoint, only : chkp_t
  use json_module, only : json_file, json_core, json_value
  use json_utils, only : json_get, json_get_or_default, json_extract_item
  use profiler, only : profiler_start_region, profiler_end_region
  use user_intf, only : user_t
  use time_step_controller, only : time_step_controller_t
  use ax_product, only : ax_t, ax_helm_factory
  use coefs, only : coef_t
  use euler_residual, only : euler_rhs_t, euler_rhs_factory
  use neko_config, only : neko_bcknd_device
  use runge_kutta_time_scheme, only : runge_kutta_time_scheme_t
  use bc_list, only : bc_list_t
  use bc, only : bc_t
  use utils, only : neko_error, neko_type_error
  use logger, only : log_size, neko_log
  use time_state, only : time_state_t
  use compressible_ops_cpu, only : compressible_ops_cpu_update_uvw, &
       compressible_ops_cpu_update_mxyz_p_ruvw, &
       compressible_ops_cpu_update_e
  use compressible_ops_device, only : compressible_ops_device_update_uvw, &
       compressible_ops_device_update_mxyz_p_ruvw, &
       compressible_ops_device_update_e
  use mpi_f08, only : mpi_allreduce, mpi_integer, mpi_max
  use regularization, only : regularization_t, regularization_factory
  implicit none
  private
 
  type, public, extends(fluid_scheme_compressible_t) &
       :: fluid_scheme_compressible_euler_t
     type(field_t) :: rho_res, m_x_res, m_y_res, m_z_res, m_e_res
     type(field_t) :: drho, dm_x, dm_y, dm_z, de
     type(field_t) :: h
     real(kind=rp) :: c_avisc_low
     class(advection_t), allocatable :: adv
     class(ax_t), allocatable :: ax
     class(euler_rhs_t), allocatable :: euler_rhs
     type(runge_kutta_time_scheme_t) :: rk_scheme
 
     class(regularization_t), allocatable :: regularization
 
     ! List of boundary conditions for velocity
     type(bc_list_t) :: bcs_density
   contains
     procedure, pass(this) :: init => fluid_scheme_compressible_euler_init
     procedure, pass(this) :: free => fluid_scheme_compressible_euler_free
     procedure, pass(this) :: step => fluid_scheme_compressible_euler_step
     procedure, pass(this) :: restart => fluid_scheme_compressible_euler_restart
     procedure, pass(this) :: setup_bcs &
          => fluid_scheme_compressible_euler_setup_bcs
     procedure, pass(this) :: compute_h
     procedure, pass(this), private :: setup_regularization
  end type fluid_scheme_compressible_euler_t
 
  interface
 
     module subroutine density_bc_factory(object, scheme, json, coef, user)
       class(bc_t), pointer, intent(inout) :: object
       type(fluid_scheme_compressible_euler_t), intent(in) :: scheme
       type(json_file), intent(inout) :: json
       type(coef_t), intent(in) :: coef
       type(user_t), intent(in) :: user
     end subroutine density_bc_factory
  end interface
 
  interface
 
     module subroutine pressure_bc_factory(object, scheme, json, coef, user)
       class(bc_t), pointer, intent(inout) :: object
       type(fluid_scheme_compressible_euler_t), intent(inout) :: scheme
       type(json_file), intent(inout) :: json
       type(coef_t), intent(in) :: coef
       type(user_t), intent(in) :: user
     end subroutine pressure_bc_factory
  end interface
 
  interface
 
     module subroutine velocity_bc_factory(object, scheme, json, coef, user)
       class(bc_t), pointer, intent(inout) :: object
       type(fluid_scheme_compressible_euler_t), intent(in) :: scheme
       type(json_file), intent(inout) :: json
       type(coef_t), intent(in) :: coef
       type(user_t), intent(in) :: user
     end subroutine velocity_bc_factory
  end interface
 
contains
  subroutine fluid_scheme_compressible_euler_init(this, msh, lx, params, user, &
       chkp)
    class(fluid_scheme_compressible_euler_t), target, intent(inout) :: this
    type(mesh_t), target, intent(inout) :: msh
    integer, intent(in) :: lx
    type(json_file), target, intent(inout) :: params
    type(user_t), target, intent(in) :: user
    type(chkp_t), target, intent(inout) :: chkp
    character(len=12), parameter :: scheme = 'compressible'
    integer :: rk_order
 
    call this%free()
 
    ! Initialize base class
    call this%scheme_init(msh, lx, params, scheme, user)
 
    call euler_rhs_factory(this%euler_rhs)
 
    associate(xh_lx => this%Xh%lx, xh_ly => this%Xh%ly, xh_lz => this%Xh%lz, &
         dm_xh => this%dm_Xh, nelv => this%msh%nelv)
 
      call this%drho%init(dm_xh, 'drho')
      call this%dm_x%init(dm_xh, 'dm_x')
      call this%dm_y%init(dm_xh, 'dm_y')
      call this%dm_z%init(dm_xh, 'dm_z')
      call this%dE%init(dm_xh, 'dE')
      call this%h%init(dm_xh, 'h')
 
    end associate
 
    if (neko_bcknd_device .eq. 1) then
       associate(p => this%p, rho => this%rho, &
            u => this%u, v => this%v, w => this%w, &
            m_x => this%m_x, m_y => this%m_y, m_z => this%m_z, &
            effective_visc => this%effective_visc)
         call device_memcpy(p%x, p%x_d, p%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(rho%x, rho%x_d, rho%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(u%x, u%x_d, u%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(v%x, v%x_d, v%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(w%x, w%x_d, w%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(m_x%x, m_x%x_d, m_x%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(m_y%x, m_y%x_d, m_y%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(m_z%x, m_z%x_d, m_z%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(effective_visc%x, effective_visc%x_d, &
              effective_visc%dof%size(), host_to_device, sync = .false.)
       end associate
    end if
 
    ! Initialize the diffusion operator
    call ax_helm_factory(this%Ax, full_formulation = .false.)
 
    ! Compute h
    call this%compute_h()
 
    ! Initialize regularization
    call this%setup_regularization(params)
 
    ! Initialize Runge-Kutta scheme
    call json_get_or_default(params, 'case.numerics.time_order', rk_order, 4)
    call this%rk_scheme%init(rk_order)
 
    call neko_log%section("Fluid boundary conditions")
    ! Set up boundary conditions
    call this%setup_bcs(user, params)
    call neko_log%end_section()
 
  end subroutine fluid_scheme_compressible_euler_init
 
  subroutine fluid_scheme_compressible_euler_free(this)
    class(fluid_scheme_compressible_euler_t), intent(inout) :: this
 
    call this%scheme_free()
 
    if (allocated(this%Ax)) then
       deallocate(this%Ax)
    end if
 
    if (allocated(this%euler_rhs)) then
       deallocate(this%euler_rhs)
    end if
 
    call this%drho%free()
    call this%dm_x%free()
    call this%dm_y%free()
    call this%dm_z%free()
    call this%dE%free()
    call this%h%free()
 
    if (allocated(this%regularization)) then
       call this%regularization%free()
       deallocate(this%regularization)
    end if
 
    call this%bcs_density%free()
 
  end subroutine fluid_scheme_compressible_euler_free
 
  subroutine fluid_scheme_compressible_euler_step(this, time, dt_controller)
    use entropy_viscosity, only : entropy_viscosity_t
    class(fluid_scheme_compressible_euler_t), target, intent(inout) :: this
    type(time_state_t), intent(in) :: time
    type(time_step_controller_t), intent(in) :: dt_controller
    type(field_t), pointer :: temp
    integer :: temp_indices(1)
    ! number of degrees of freedom
    integer :: n
    integer :: i
    class(bc_t), pointer :: b
 
    n = this%dm_Xh%size()
    call neko_scratch_registry%request_field(temp, temp_indices(1), .false.)
    b => null()
 
    call profiler_start_region('Fluid compressible', 1)
    associate(u => this%u, v => this%v, w => this%w, p => this%p, &
         m_x=> this%m_x, m_y => this%m_y, m_z => this%m_z, &
         xh => this%Xh, msh => this%msh, ax => this%Ax, &
         c_xh => this%c_Xh, dm_xh => this%dm_Xh, gs_xh => this%gs_Xh, &
         e => this%E, rho => this%rho, mu => this%mu, &
         f_x => this%f_x, f_y => this%f_y, f_z => this%f_z, &
         drho => this%drho, dm_x => this%dm_x, dm_y => this%dm_y, &
         dm_z => this%dm_z, de => this%dE, &
         euler_rhs => this%euler_rhs, h => this%h, &
         t => time%t, tstep => time%tstep, dt => time%dt, &
         c_avisc_low => this%c_avisc_low, rk_scheme => this%rk_scheme)
 
      ! Compute artificial viscosity
      call this%regularization%compute(time, time%tstep, time%dt)
 
      ! Execute RHS step with effective viscosity field
      call euler_rhs%step(rho, m_x, m_y, m_z, e, &
           p, u, v, w, ax, &
           c_xh, gs_xh, h, this%effective_visc, &
           rk_scheme, dt)
 
      call this%bcs_density%apply(rho, time)
 
      ! Update u, v, w
      if (neko_bcknd_device .eq. 1) then
         call compressible_ops_device_update_uvw(u%x_d, v%x_d, w%x_d, &
              m_x%x_d, m_y%x_d, m_z%x_d, rho%x_d, n)
      else
         call compressible_ops_cpu_update_uvw(u%x, v%x, w%x, &
              m_x%x, m_y%x, m_z%x, rho%x, n)
      end if
 
      call this%bcs_vel%apply_vector(u%x, v%x, w%x, &
           dm_xh%size(), time, strong = .true.)
 
      if (neko_bcknd_device .eq. 1) then
         call compressible_ops_device_update_mxyz_p_ruvw(m_x%x_d, m_y%x_d, &
              m_z%x_d, p%x_d, temp%x_d, u%x_d, v%x_d, w%x_d, e%x_d, &
              rho%x_d, this%gamma, n)
      else
         call compressible_ops_cpu_update_mxyz_p_ruvw(m_x%x, m_y%x, m_z%x, &
              p%x, temp%x, u%x, v%x, w%x, e%x, rho%x, this%gamma, n)
      end if
 
      call this%bcs_prs%apply(p, time)
 
 
      if (neko_bcknd_device .eq. 1) then
         call compressible_ops_device_update_e(e%x_d, p%x_d, &
              temp%x_d, this%gamma, n)
      else
         call compressible_ops_cpu_update_e(e%x, p%x, temp%x, this%gamma, n)
      end if
 
 
      call this%compute_entropy()
 
      if (allocated(this%regularization)) then
         select type (reg => this%regularization)
         type is (entropy_viscosity_t)
            call reg%update_lag()
         end select
      end if
 
      call this%compute_max_wave_speed()
 
      do i = 1, this%bcs_vel%size()
         b => this%bcs_vel%get(i)
         b%updated = .false.
      end do
 
      do i = 1, this%bcs_prs%size()
         b => this%bcs_prs%get(i)
         b%updated = .false.
      end do
 
      do i = 1, this%bcs_density%size()
         b => this%bcs_density%get(i)
         b%updated = .false.
      end do
      nullify(b)
 
    end associate
    call profiler_end_region('Fluid compressible', 1)
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine fluid_scheme_compressible_euler_step
 
  subroutine fluid_scheme_compressible_euler_setup_bcs(this, user, params)
    class(fluid_scheme_compressible_euler_t), target, intent(inout) :: this
    type(user_t), target, intent(in) :: user
    type(json_file), intent(inout) :: params
    integer :: i, n_bcs, zone_index, j, zone_size, global_zone_size, ierr
    class(bc_t), pointer :: bc_i
    type(json_core) :: core
    type(json_value), pointer :: bc_object
    type(json_file) :: bc_subdict
    logical :: found
    integer, allocatable :: zone_indices(:)
    character(len=LOG_SIZE) :: log_buf
 
    ! Process boundary conditions
    if (params%valid_path('case.fluid.boundary_conditions')) then
       call params%info('case.fluid.boundary_conditions', n_children = n_bcs)
       call params%get_core(core)
       call params%get('case.fluid.boundary_conditions', bc_object, found)
 
       !
       ! Velocity bcs
       !
       call this%bcs_vel%init(n_bcs)
 
       do i = 1, n_bcs
          ! Extract BC configuration
          call json_extract_item(core, bc_object, i, bc_subdict)
          call json_get(bc_subdict, "zone_indices", zone_indices)
 
          ! Validate zones
          do j = 1, size(zone_indices)
             zone_size = this%msh%labeled_zones(zone_indices(j))%size
             call mpi_allreduce(zone_size, global_zone_size, 1, &
                  mpi_integer, mpi_max, neko_comm, ierr)
 
             if (global_zone_size .eq. 0) then
                write(log_buf, '(A,I0,A)') "Error: Zone ", zone_indices(j), &
                     " has zero size"
                call neko_error(log_buf)
             end if
          end do
 
          ! Create BC
          bc_i => null()
          call velocity_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
 
          ! Add to appropriate lists
          if (associated(bc_i)) then
             call this%bcs_vel%append(bc_i)
          end if
       end do
 
       !
       ! Pressure bcs
       !
       call this%bcs_prs%init(n_bcs)
 
       do i = 1, n_bcs
          ! Create a new json containing just the subdict for this bc
          call json_extract_item(core, bc_object, i, bc_subdict)
          bc_i => null()
          call pressure_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
 
          ! Not all bcs require an allocation for pressure in particular,
          ! so we check.
          if (associated(bc_i)) then
             call this%bcs_prs%append(bc_i)
          end if
       end do
 
       !
       ! Density bcs
       !
       call this%bcs_density%init(n_bcs)
 
       do i = 1, n_bcs
          ! Create a new json containing just the subdict for this bc
          call json_extract_item(core, bc_object, i, bc_subdict)
          bc_i => null()
          call density_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
 
          ! Not all bcs require an allocation for pressure in particular,
          ! so we check.
          if (associated(bc_i)) then
             call this%bcs_density%append(bc_i)
          end if
       end do
    else
       ! Check that there are no labeled zones, i.e. all are periodic.
       do i = 1, size(this%msh%labeled_zones)
          if (this%msh%labeled_zones(i)%size .gt. 0) then
             call neko_error("No boundary_conditions entry in the case file!")
          end if
       end do
 
       ! For a pure periodic case, we still need to initilise the bc lists
       ! to a zero size to avoid issues with apply() in step()
       call this%bcs_prs%init()
       call this%bcs_vel%init()
       call this%bcs_density%init()
 
    end if
  end subroutine fluid_scheme_compressible_euler_setup_bcs
 
  subroutine compute_h(this)
    class(fluid_scheme_compressible_euler_t), intent(inout) :: this
    integer :: lx, ly, lz
 
    lx = this%c_Xh%Xh%lx
    ly = this%c_Xh%Xh%ly
    lz = this%c_Xh%Xh%lz
    call compute_h_cpu(this%h%x, this%c_Xh%dof%x, this%c_Xh%dof%y, &
         this%c_Xh%dof%z, lx, ly, lz, this%c_Xh%msh%nelv)
 
    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_Xh%op(this%h, gs_op_add)
       call device_col2(this%h%x_d, this%c_Xh%mult_d, this%h%dof%size())
    else
       call this%gs_Xh%op(this%h, gs_op_add)
       call col2(this%h%x, this%c_Xh%mult, this%h%dof%size())
    end if
 
  end subroutine compute_h
 
  subroutine compute_h_cpu(h, x, y, z, lx, ly, lz, nelv)
    integer, intent(in) :: lx, ly, lz, nelv
    real(kind=rp), intent(out) :: h(lx, ly, lz, nelv)
    real(kind=rp), intent(in) :: x(lx, ly, lz, nelv)
    real(kind=rp), intent(in) :: y(lx, ly, lz, nelv)
    real(kind=rp), intent(in) :: z(lx, ly, lz, nelv)
    integer :: e, i, j, k
    integer :: im, ip, jm, jp, km, kp
    real(kind=rp) :: di, dj, dk
 
    !$omp parallel do private(i, j, k, im, ip, jm, jp, km, kp, di, dj, dk)
    do e = 1, nelv
       do k = 1, lz
          km = max(1, k - 1)
          kp = min(lz, k + 1)
          do j = 1, ly
             jm = max(1, j - 1)
             jp = min(ly, j + 1)
             do i = 1, lx
                im = max(1, i - 1)
                ip = min(lx, i + 1)
 
                di = (x(ip, j, k, e) - x(im, j, k, e))**2 &
                     + (y(ip, j, k, e) - y(im, j, k, e))**2 &
                     + (z(ip, j, k, e) - z(im, j, k, e))**2
 
                dj = (x(i, jp, k, e) - x(i, jm, k, e))**2 &
                     + (y(i, jp, k, e) - y(i, jm, k, e))**2 &
                     + (z(i, jp, k, e) - z(i, jm, k, e))**2
 
                dk = (x(i, j, kp, e) - x(i, j, km, e))**2 &
                     + (y(i, j, kp, e) - y(i, j, km, e))**2 &
                     + (z(i, j, kp, e) - z(i, j, km, e))**2
 
                di = sqrt(di) / (ip - im)
                dj = sqrt(dj) / (jp - jm)
                dk = sqrt(dk) / (kp - km)
                h(i,j,k,e) = (di * dj * dk)**(1.0_rp / 3.0_rp)
             end do
          end do
       end do
    end do
    !$omp end parallel do
  end subroutine compute_h_cpu
 
  subroutine fluid_scheme_compressible_euler_restart(this, chkp)
    class(fluid_scheme_compressible_euler_t), target, intent(inout) :: this
    type(chkp_t), intent(inout) :: chkp
  end subroutine fluid_scheme_compressible_euler_restart
 
  subroutine setup_regularization(this, params)
    use entropy_viscosity, only : entropy_viscosity_t, &
         entropy_viscosity_set_fields
    class(fluid_scheme_compressible_euler_t), target, intent(inout) :: this
    type(json_file), intent(inout) :: params
    type(json_file) :: reg_json
    type(json_core) :: json_core_inst
    type(json_value), pointer :: reg_params
    character(len=:), allocatable :: buffer
    real(kind=rp) :: c_avisc_entropy_val
    character(len=:), allocatable :: regularization_type
 
    call json_get_or_default(params, 'case.numerics.c_avisc_low', &
         this%c_avisc_low, 0.5_rp)
    call json_get_or_default(params, 'case.numerics.c_avisc_entropy', &
         c_avisc_entropy_val, 1.0_rp)
 
    call json_core_inst%initialize()
    call json_core_inst%create_object(reg_params, '')
    call json_core_inst%add(reg_params, 'c_avisc_entropy', c_avisc_entropy_val)
    call json_core_inst%add(reg_params, 'c_avisc_low', this%c_avisc_low)
    call json_core_inst%print_to_string(reg_params, buffer)
    call json_core_inst%destroy(reg_params)
 
    call reg_json%initialize()
    call reg_json%load_from_string(buffer)
 
    regularization_type = 'entropy_viscosity'
 
    call regularization_factory(this%regularization, regularization_type, &
         reg_json, this%c_Xh, this%dm_Xh, this%effective_visc)
 
    select type (reg => this%regularization)
    type is (entropy_viscosity_t)
       call entropy_viscosity_set_fields(reg, this%S, this%u, this%v, this%w, &
            this%h, this%max_wave_speed, this%msh, this%Xh, this%gs_Xh)
    end select
 
    call reg_json%destroy()
 
  end subroutine setup_regularization
 
end module fluid_scheme_compressible_euler