fluid_scheme_compressible_ns.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_ns
  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 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_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_allocator
  use coefs, only : coef_t
  use compressible_residual, only : compressible_rhs_t, compressible_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, &
       compressible_ops_device_update_temperature
  use neko_config, only : neko_bcknd_device
  use mpi_f08, only : mpi_allreduce, mpi_integer, mpi_max
  use viscous_regularization, only : viscous_regularization_t, &
       viscous_regularization_factory
  use vector_bc_projector, only : vector_bc_projector_t, &
       coupled_vector_bc_projector_t
  implicit none
  private
 
  type, public, extends(fluid_scheme_compressible_t) &
       :: fluid_scheme_compressible_ns_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
     class(advection_t), allocatable :: adv
     class(ax_t), allocatable :: ax
     class(ax_t), allocatable :: ax_stress
     class(compressible_rhs_t), allocatable :: compressible_rhs
     type(runge_kutta_time_scheme_t) :: rk_scheme
 
     class(viscous_regularization_t), allocatable :: viscous_regularization
     logical :: if_viscous_regularization
 
     type(coupled_vector_bc_projector_t):: bcs_vel_projector
     type(bc_list_t) :: bcs_density
 
   contains
     procedure, pass(this) :: init => fluid_scheme_compressible_ns_init
     procedure, pass(this) :: free => fluid_scheme_compressible_ns_free
     procedure, pass(this) :: step => fluid_scheme_compressible_ns_step
     procedure, pass(this) :: restart => fluid_scheme_compressible_ns_restart
     procedure, pass(this) :: setup_bcs &
          => fluid_scheme_compressible_ns_setup_bcs
     procedure, pass(this), private :: setup_viscous_regularization
  end type fluid_scheme_compressible_ns_t
 
  interface
 
     module subroutine density_bc_factory(object, scheme, json, coef, user)
       class(bc_t), pointer, intent(inout) :: object
       type(fluid_scheme_compressible_ns_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_ns_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_ns_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_ns_init(this, msh, lx, params, user, &
       chkp)
    class(fluid_scheme_compressible_ns_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 compressible_rhs_factory(this%compressible_rhs, this%gamma)
 
    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')
 
    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, &
            artificial_visc => this%artificial_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(artificial_visc%x, artificial_visc%x_d, &
              artificial_visc%dof%size(), host_to_device, sync = .false.)
       end associate
    end if
 
    ! Initialize the diffusion operators
    if (this%c_Xh%geo_compression) then
       call ax_helm_allocator(this%Ax, type_name = "standard_compr")
    else
       call ax_helm_allocator(this%Ax, type_name = "standard")
    end if
    call ax_helm_allocator(this%Ax_stress, type_name = "full")
 
    ! Initialize the velocity BC projector. The coupled projector builds the
    ! local bases required by non-axis aligned mixed velocity conditions.
    call this%bcs_vel_projector%init(this%c_Xh)
 
    ! Initialize viscous regularization
    call this%setup_viscous_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_ns_init
 
  subroutine fluid_scheme_compressible_ns_free(this)
    class(fluid_scheme_compressible_ns_t), intent(inout) :: this
    class(bc_t), pointer :: bc
    integer :: i
 
    if (allocated(this%Ax)) then
       call this%Ax%free()
       deallocate(this%Ax)
    end if
 
    if (allocated(this%Ax_stress)) then
       call this%Ax_stress%free()
       deallocate(this%Ax_stress)
    end if
 
    call this%scheme_free()
 
    if (allocated(this%compressible_rhs)) then
       deallocate(this%compressible_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%rk_scheme%free()
 
    if (allocated(this%viscous_regularization)) then
       call this%viscous_regularization%free()
       deallocate(this%viscous_regularization)
    end if
 
    do i = 1, this%bcs_density%size()
       bc => this%bcs_density%get(i)
       if (associated(bc)) then
          call bc%free()
          deallocate(bc)
       end if
    end do
    call this%bcs_density%free()
    call this%bcs_vel_projector%free()
 
  end subroutine fluid_scheme_compressible_ns_free
 
  subroutine fluid_scheme_compressible_ns_step(this, time, dt_controller)
    class(fluid_scheme_compressible_ns_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, &
         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, &
         compressible_rhs => this%compressible_rhs, &
         t => time%t, tstep => time%tstep, dt => time%dt, &
         rk_scheme => this%rk_scheme)
 
      if (allocated(this%viscous_regularization)) then
         call this%viscous_regularization%update(this%artificial_visc)
      end if
 
      ! Refresh user-specified physical viscosity/conductivity before RHS.
      call this%update_material_properties(time)
 
      ! Execute RHS step with artificial viscosity field
      call compressible_rhs%step(rho, m_x, m_y, m_z, e, &
           p, u, v, w, this%Ax, &
           this%Ax_stress, c_xh, gs_xh, this%artificial_visc, this%mu, &
           this%kappa, this%bcs_vel, time, rk_scheme, real(dt, kind=rp))
 
      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
 
      if (neko_bcknd_device .eq. 1) then
         call compressible_ops_device_update_temperature( &
              this%temperature%x_d, p%x_d, rho%x_d, this%gamma, n)
      else
         !OCL NORECURRENCE, NOVREC, NOALIAS
         !DIR$ CONCURRENT
         !DIR$ IVDEP
         !GCC$ ivdep
         !$omp parallel do simd
         do i = 1, n
            this%temperature%x(i,1,1,1) = p%x(i,1,1,1) / &
                 (rho%x(i,1,1,1) * (this%gamma - 1.0_rp))
         end do
         !$omp end parallel do simd
      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_ns_step
 
  subroutine fluid_scheme_compressible_ns_setup_bcs(this, user, params)
    class(fluid_scheme_compressible_ns_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_projector%mark(bc_i)
             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
 
    call this%bcs_vel_projector%finalize(rebuild_mask = .true.)
  end subroutine fluid_scheme_compressible_ns_setup_bcs
 
  subroutine fluid_scheme_compressible_ns_restart(this, chkp)
    class(fluid_scheme_compressible_ns_t), target, intent(inout) :: this
    type(chkp_t), intent(inout) :: chkp
  end subroutine fluid_scheme_compressible_ns_restart
 
  subroutine setup_viscous_regularization(this, params)
    class(fluid_scheme_compressible_ns_t), target, intent(inout) :: this
    type(json_file), intent(inout) :: params
    type(json_file) :: reg_json
    logical :: found
    character(len=:), allocatable :: viscous_regularization_type
 
    found = .false.
    if (this%params%valid_path('case.fluid.viscous_regularization')) then
       call json_get(params, 'case.fluid.viscous_regularization', &
            reg_json)
       found = .true.
    end if
 
    if (.not. found) then
       ! No viscous_regularization specified, so we skip setup
       this%if_viscous_regularization = .false.
       return
    else
       this%if_viscous_regularization = .true.
    end if
 
    call json_get(reg_json, 'type', viscous_regularization_type)
    call viscous_regularization_factory(this%viscous_regularization, &
         viscous_regularization_type, &
         reg_json, this%c_Xh, this%dm_Xh)
 
    call reg_json%destroy()
 
    if (allocated(viscous_regularization_type)) &
         deallocate(viscous_regularization_type)
 
  end subroutine setup_viscous_regularization
 
end module fluid_scheme_compressible_ns