d8/d6d/euler__res__cpu_8f90_source.html

! 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.

!

!? It handles the time advancement of primitive variables using Runge-Kutta methods

!? and evaluates the right-hand side terms of the Euler equations including artificial viscosity.

module euler_res_cpu

  use euler_residual, only : euler_rhs_t

  use field, only : field_t

  use ax_product, only : ax_t

  use coefs, only : coef_t

  use gather_scatter, only : gs_t

  use num_types, only : rp

  use operators, only : div, rotate_cyc

  use gs_ops, only : gs_op_add

  use scratch_registry, only : neko_scratch_registry

  use runge_kutta_time_scheme, only : runge_kutta_time_scheme_t

  use field_list, only : field_list_t

  implicit none

  private


  type, public, extends(euler_rhs_t) :: euler_res_cpu_t

   contains

     procedure, nopass :: step => advance_primitive_variables_cpu

     procedure, nopass :: evaluate_rhs_cpu

  end type euler_res_cpu_t


contains


  subroutine advance_primitive_variables_cpu(rho_field, m_x, m_y, m_z, &

       E, p, u, v, w, Ax, &

       coef, gs, h, effective_visc, rk_scheme, dt)

    type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E

    type(field_t), intent(in) :: p, u, v, w, h, effective_visc

    class(ax_t), intent(inout) :: Ax

    type(coef_t), intent(inout) :: coef

    type(gs_t), intent(inout) :: gs

    class(runge_kutta_time_scheme_t), intent(in) :: rk_scheme

    real(kind=rp), intent(in) :: dt

    integer :: n, s, i, j, k

    type(field_t), pointer :: k_rho_1, k_rho_2, k_rho_3, k_rho_4, &

         k_m_x_1, k_m_x_2, k_m_x_3, k_m_x_4, &

         k_m_y_1, k_m_y_2, k_m_y_3, k_m_y_4, &

         k_m_z_1, k_m_z_2, k_m_z_3, k_m_z_4, &

         k_E_1, k_E_2, k_E_3, k_E_4, &

         temp_rho, temp_m_x, temp_m_y, temp_m_z, temp_E

    integer :: tmp_indices(25)

    type(field_list_t) :: k_rho, k_m_x, k_m_y, k_m_z, k_E

    ! These contiguous pointers are necessary to ensure that

    ! the Fujitsu compiler properly vectorizes the loop nests

    real(kind=rp), contiguous, pointer :: k_rho_ptr(:,:,:,:)

    real(kind=rp), contiguous, pointer :: k_m_x_ptr(:,:,:,:)

    real(kind=rp), contiguous, pointer :: k_m_y_ptr(:,:,:,:)

    real(kind=rp), contiguous, pointer :: k_m_z_ptr(:,:,:,:)

    real(kind=rp), contiguous, pointer :: k_e_ptr(:,:,:,:)


    n = p%dof%size()

    s = rk_scheme%order

    call neko_scratch_registry%request_field(k_rho_1, tmp_indices(1), .true.)

    call neko_scratch_registry%request_field(k_rho_2, tmp_indices(2), .true.)

    call neko_scratch_registry%request_field(k_rho_3, tmp_indices(3), .true.)

    call neko_scratch_registry%request_field(k_rho_4, tmp_indices(4), .true.)

    call neko_scratch_registry%request_field(k_m_x_1, tmp_indices(5), .true.)

    call neko_scratch_registry%request_field(k_m_x_2, tmp_indices(6), .true.)

    call neko_scratch_registry%request_field(k_m_x_3, tmp_indices(7), .true.)

    call neko_scratch_registry%request_field(k_m_x_4, tmp_indices(8), .true.)

    call neko_scratch_registry%request_field(k_m_y_1, tmp_indices(9), .true.)

    call neko_scratch_registry%request_field(k_m_y_2, tmp_indices(10), .true.)

    call neko_scratch_registry%request_field(k_m_y_3, tmp_indices(11), .true.)

    call neko_scratch_registry%request_field(k_m_y_4, tmp_indices(12), .true.)

    call neko_scratch_registry%request_field(k_m_z_1, tmp_indices(13), .true.)

    call neko_scratch_registry%request_field(k_m_z_2, tmp_indices(14), .true.)

    call neko_scratch_registry%request_field(k_m_z_3, tmp_indices(15), .true.)

    call neko_scratch_registry%request_field(k_m_z_4, tmp_indices(16), .true.)

    call neko_scratch_registry%request_field(k_e_1, tmp_indices(17), .true.)

    call neko_scratch_registry%request_field(k_e_2, tmp_indices(18), .true.)

    call neko_scratch_registry%request_field(k_e_3, tmp_indices(19), .true.)

    call neko_scratch_registry%request_field(k_e_4, tmp_indices(20), .true.)

    call neko_scratch_registry%request_field(temp_rho, tmp_indices(21), .false.)

    call neko_scratch_registry%request_field(temp_m_x, tmp_indices(22), .false.)

    call neko_scratch_registry%request_field(temp_m_y, tmp_indices(23), .false.)

    call neko_scratch_registry%request_field(temp_m_z, tmp_indices(24), .false.)

    call neko_scratch_registry%request_field(temp_e, tmp_indices(25), .false.)


    ! Initialize Runge-Kutta stage variables for each conserved quantity

    call k_rho%init(4)

    call k_rho%assign(1, k_rho_1)

    call k_rho%assign(2, k_rho_2)

    call k_rho%assign(3, k_rho_3)

    call k_rho%assign(4, k_rho_4)

    call k_m_x%init(4)

    call k_m_x%assign(1, k_m_x_1)

    call k_m_x%assign(2, k_m_x_2)

    call k_m_x%assign(3, k_m_x_3)

    call k_m_x%assign(4, k_m_x_4)

    call k_m_y%init(4)

    call k_m_y%assign(1, k_m_y_1)

    call k_m_y%assign(2, k_m_y_2)

    call k_m_y%assign(3, k_m_y_3)

    call k_m_y%assign(4, k_m_y_4)

    call k_m_z%init(4)

    call k_m_z%assign(1, k_m_z_1)

    call k_m_z%assign(2, k_m_z_2)

    call k_m_z%assign(3, k_m_z_3)

    call k_m_z%assign(4, k_m_z_4)

    call k_e%init(4)

    call k_e%assign(1, k_e_1)

    call k_e%assign(2, k_e_2)

    call k_e%assign(3, k_e_3)

    call k_e%assign(4, k_e_4)


    ! Loop over Runge-Kutta stages

    do i = 1, s

       ! Copy current solution state to temporary arrays for this RK stage

       do concurrent(k = 1:n)

          temp_rho%x(k,1,1,1) = rho_field%x(k,1,1,1)

          temp_m_x%x(k,1,1,1) = m_x%x(k,1,1,1)

          temp_m_y%x(k,1,1,1) = m_y%x(k,1,1,1)

          temp_m_z%x(k,1,1,1) = m_z%x(k,1,1,1)

          temp_e%x(k,1,1,1) = e%x(k,1,1,1)

       end do


       ! Accumulate previous stage contributions using RK coefficients

       do j = 1, i-1

          k_rho_ptr => k_rho%items(j)%ptr%x

          k_m_x_ptr => k_m_x%items(j)%ptr%x

          k_m_y_ptr => k_m_y%items(j)%ptr%x

          k_m_z_ptr => k_m_z%items(j)%ptr%x

          k_e_ptr => k_e%items(j)%ptr%x

          do concurrent(k = 1:n)

             temp_rho%x(k,1,1,1) = temp_rho%x(k,1,1,1) &

                  + dt * rk_scheme%coeffs_A(i, j) * k_rho_ptr(k,1,1,1)

             temp_m_x%x(k,1,1,1) = temp_m_x%x(k,1,1,1) &

                  + dt * rk_scheme%coeffs_A(i, j) * k_m_x_ptr(k,1,1,1)

             temp_m_y%x(k,1,1,1) = temp_m_y%x(k,1,1,1) &

                  + dt * rk_scheme%coeffs_A(i, j) * k_m_y_ptr(k,1,1,1)

             temp_m_z%x(k,1,1,1) = temp_m_z%x(k,1,1,1) &

                  + dt * rk_scheme%coeffs_A(i, j) * k_m_z_ptr(k,1,1,1)

             temp_e%x(k,1,1,1) = temp_e%x(k,1,1,1) &

                  + dt * rk_scheme%coeffs_A(i, j) * k_e_ptr(k,1,1,1)

          end do

       end do


       ! Evaluate RHS terms for current stage using intermediate solution values

       call evaluate_rhs_cpu(k_rho%items(i)%ptr, k_m_x%items(i)%ptr, &

            k_m_y%items(i)%ptr, k_m_z%items(i)%ptr, &

            k_e%items(i)%ptr, &

            temp_rho, temp_m_x, temp_m_y, temp_m_z, temp_e, &

            p, u, v, w, ax, &

            coef, gs, h, effective_visc)

    end do


    ! Update the solution

    do i = 1, s

       k_rho_ptr => k_rho%items(i)%ptr%x

       k_m_x_ptr => k_m_x%items(i)%ptr%x

       k_m_y_ptr => k_m_y%items(i)%ptr%x

       k_m_z_ptr => k_m_z%items(i)%ptr%x

       k_e_ptr => k_e%items(i)%ptr%x

       do concurrent(k = 1:n)

          rho_field%x(k,1,1,1) = rho_field%x(k,1,1,1) &

               + dt * rk_scheme%coeffs_b(i) * k_rho_ptr(k,1,1,1)

          m_x%x(k,1,1,1) = m_x%x(k,1,1,1) &

               + dt * rk_scheme%coeffs_b(i) * k_m_x_ptr(k,1,1,1)

          m_y%x(k,1,1,1) = m_y%x(k,1,1,1) &

               + dt * rk_scheme%coeffs_b(i) * k_m_y_ptr(k,1,1,1)

          m_z%x(k,1,1,1) = m_z%x(k,1,1,1) &

               + dt * rk_scheme%coeffs_b(i) * k_m_z_ptr(k,1,1,1)

          e%x(k,1,1,1) = e%x(k,1,1,1) &

               + dt * rk_scheme%coeffs_b(i) * k_e_ptr(k,1,1,1)

       end do

    end do


    call neko_scratch_registry%relinquish_field(tmp_indices)


  end subroutine advance_primitive_variables_cpu


  subroutine evaluate_rhs_cpu(rhs_rho_field, rhs_m_x, rhs_m_y, rhs_m_z, rhs_E, &

       rho_field, m_x, m_y, m_z, E, p, u, v, w, Ax, &

       coef, gs, h, effective_visc)

    type(field_t), intent(inout) :: rhs_rho_field, &

         rhs_m_x, rhs_m_y, rhs_m_z, rhs_e

    type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E

    type(field_t), intent(in) :: p, u, v, w, h, effective_visc

    class(ax_t), intent(inout) :: Ax

    type(coef_t), intent(inout) :: coef

    type(gs_t), intent(inout) :: gs

    integer :: i, n

    type(field_t), pointer :: f_x, f_y, f_z, &

         visc_rho, visc_m_x, visc_m_y, visc_m_z, visc_E

    integer :: tmp_indices(8)


    n = coef%dof%size()

    call neko_scratch_registry%request_field(f_x, tmp_indices(1), .false.)

    call neko_scratch_registry%request_field(f_y, tmp_indices(2), .false.)

    call neko_scratch_registry%request_field(f_z, tmp_indices(3), .false.)


    ! Compute density flux divergence

    call div(rhs_rho_field%x, m_x%x, m_y%x, m_z%x, coef)


    ! Compute momentum flux divergences

    ! m_x

    do concurrent(i = 1:n)

       f_x%x(i,1,1,1) = m_x%x(i,1,1,1) * m_x%x(i,1,1,1) / &

            rho_field%x(i, 1, 1, 1) + p%x(i,1,1,1)

       f_y%x(i,1,1,1) = m_x%x(i,1,1,1) * m_y%x(i,1,1,1) / &

            rho_field%x(i, 1, 1, 1)

       f_z%x(i,1,1,1) = m_x%x(i,1,1,1) * m_z%x(i,1,1,1) / &

            rho_field%x(i, 1, 1, 1)

    end do

    call div(rhs_m_x%x, f_x%x, f_y%x, f_z%x, coef)

    ! m_y

    do concurrent(i = 1:n)

       f_x%x(i,1,1,1) = m_y%x(i,1,1,1) * m_x%x(i,1,1,1) / &

            rho_field%x(i, 1, 1, 1)

       f_y%x(i,1,1,1) = m_y%x(i,1,1,1) * m_y%x(i,1,1,1) / &

            rho_field%x(i, 1, 1, 1) + p%x(i,1,1,1)

       f_z%x(i,1,1,1) = m_y%x(i,1,1,1) * m_z%x(i,1,1,1) / &

            rho_field%x(i, 1, 1, 1)

    end do

    call div(rhs_m_y%x, f_x%x, f_y%x, f_z%x, coef)

    ! m_z

    do concurrent(i = 1:n)

       f_x%x(i,1,1,1) = m_z%x(i,1,1,1) * m_x%x(i,1,1,1) / &

            rho_field%x(i, 1, 1, 1)

       f_y%x(i,1,1,1) = m_z%x(i,1,1,1) * m_y%x(i,1,1,1) / &

            rho_field%x(i, 1, 1, 1)

       f_z%x(i,1,1,1) = m_z%x(i,1,1,1) * m_z%x(i,1,1,1) / &

            rho_field%x(i, 1, 1, 1) + p%x(i,1,1,1)

    end do

    call div(rhs_m_z%x, f_x%x, f_y%x, f_z%x, coef)


    ! Compute energy flux divergence

    do concurrent(i = 1:n)

       f_x%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) * u%x(i,1,1,1)

       f_y%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) * v%x(i,1,1,1)

       f_z%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) * w%x(i,1,1,1)

    end do

    call div(rhs_e%x, f_x%x, f_y%x, f_z%x, coef)


    ! gs

    call gs%op(rhs_rho_field, gs_op_add)

    call rotate_cyc(rhs_m_x%x, rhs_m_y%x, rhs_m_z%x, 1, coef)

    call gs%op(rhs_m_x, gs_op_add)

    call gs%op(rhs_m_y, gs_op_add)

    call gs%op(rhs_m_z, gs_op_add)

    call rotate_cyc(rhs_m_x%x, rhs_m_y%x, rhs_m_z%x, 0, coef)

    call gs%op(rhs_e, gs_op_add)

    do concurrent(i = 1:rhs_e%dof%size())

       rhs_rho_field%x(i,1,1,1) = rhs_rho_field%x(i,1,1,1) * coef%mult(i,1,1,1)

       rhs_m_x%x(i,1,1,1) = rhs_m_x%x(i,1,1,1) * coef%mult(i,1,1,1)

       rhs_m_y%x(i,1,1,1) = rhs_m_y%x(i,1,1,1) * coef%mult(i,1,1,1)

       rhs_m_z%x(i,1,1,1) = rhs_m_z%x(i,1,1,1) * coef%mult(i,1,1,1)

       rhs_e%x(i,1,1,1) = rhs_e%x(i,1,1,1) * coef%mult(i,1,1,1)

    end do


    call neko_scratch_registry%request_field(visc_rho, tmp_indices(4), .false.)

    call neko_scratch_registry%request_field(visc_m_x, tmp_indices(5), .false.)

    call neko_scratch_registry%request_field(visc_m_y, tmp_indices(6), .false.)

    call neko_scratch_registry%request_field(visc_m_z, tmp_indices(7), .false.)

    call neko_scratch_registry%request_field(visc_e, tmp_indices(8), .false.)


    ! Set h1 coefficient to the effective viscosity for the Laplacian operator

    do concurrent(i = 1:n)

       coef%h1(i,1,1,1) = effective_visc%x(i,1,1,1)

    end do


    ! Calculate artificial diffusion with variable viscosity

    call ax%compute(visc_rho%x, rho_field%x, coef, p%msh, p%Xh)

    call ax%compute(visc_m_x%x, m_x%x, coef, p%msh, p%Xh)

    call ax%compute(visc_m_y%x, m_y%x, coef, p%msh, p%Xh)

    call ax%compute(visc_m_z%x, m_z%x, coef, p%msh, p%Xh)

    call ax%compute(visc_e%x, e%x, coef, p%msh, p%Xh)


    ! Reset h1 coefficient back to 1.0 for other operations

    do concurrent(i = 1:n)

       coef%h1(i,1,1,1) = 1.0_rp

    end do


    ! gs

    call gs%op(visc_rho, gs_op_add)

    call rotate_cyc(visc_m_x%x, visc_m_y%x, visc_m_z%x, 1, coef)

    call gs%op(visc_m_x, gs_op_add)

    call gs%op(visc_m_y, gs_op_add)

    call gs%op(visc_m_z, gs_op_add)

    call rotate_cyc(visc_m_x%x, visc_m_y%x, visc_m_z%x, 0, coef)

    call gs%op(visc_e, gs_op_add)


    ! Move div to the rhs and apply artificial viscosity

    ! The viscosity coefficient is already included in the Laplacian operator

    do concurrent(i = 1:n)

       rhs_rho_field%x(i,1,1,1) = -rhs_rho_field%x(i,1,1,1) &

            - coef%Binv(i,1,1,1) * visc_rho%x(i,1,1,1)

       rhs_m_x%x(i,1,1,1) = -rhs_m_x%x(i,1,1,1) &

            - coef%Binv(i,1,1,1) * visc_m_x%x(i,1,1,1)

       rhs_m_y%x(i,1,1,1) = -rhs_m_y%x(i,1,1,1) &

            - coef%Binv(i,1,1,1) * visc_m_y%x(i,1,1,1)

       rhs_m_z%x(i,1,1,1) = -rhs_m_z%x(i,1,1,1) &

            - coef%Binv(i,1,1,1) * visc_m_z%x(i,1,1,1)

       rhs_e%x(i,1,1,1) = -rhs_e%x(i,1,1,1) &

            - coef%Binv(i,1,1,1) * visc_e%x(i,1,1,1)

    end do


    call neko_scratch_registry%relinquish_field(tmp_indices)


  end subroutine evaluate_rhs_cpu


end module euler_res_cpu

operators::rotate_cyc
Definition operators.f90:77

ax_product
Defines a Matrix-vector product.
Definition ax.f90:34

coefs
Coefficients.
Definition coef.f90:34

euler_res_cpu
This module implements CPU-based residual calculations for the Euler equations.
Definition euler_res_cpu.f90:36

euler_res_cpu::evaluate_rhs_cpu
subroutine evaluate_rhs_cpu(rhs_rho_field, rhs_m_x, rhs_m_y, rhs_m_z, rhs_e, rho_field, m_x, m_y, m_z, e, p, u, v, w, ax, coef, gs, h, effective_visc)
Evaluates the right-hand side of the Euler equations Computes fluxes, divergence terms,...
Definition euler_res_cpu.f90:244

euler_res_cpu::advance_primitive_variables_cpu
subroutine advance_primitive_variables_cpu(rho_field, m_x, m_y, m_z, e, p, u, v, w, ax, coef, gs, h, effective_visc, rk_scheme, dt)
Advances the primitive variables (density, momentum, energy) in time using a Runge-Kutta scheme.
Definition euler_res_cpu.f90:75

euler_residual
Definition euler_res.f90:33

field_list
Definition field_list.f90:1

field
Defines a field.
Definition field.f90:34

gather_scatter
Gather-scatter.
Definition gather_scatter.f90:34

gs_ops
Defines Gather-scatter operations.
Definition gs_ops.f90:34

gs_ops::gs_op_add
integer, parameter, public gs_op_add
Definition gs_ops.f90:36

num_types
Definition num_types.f90:1

num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12

operators
Operators.
Definition operators.f90:34

operators::div
subroutine, public div(res, ux, uy, uz, coef)
Compute the divergence of a vector field.
Definition operators.f90:117

runge_kutta_time_scheme
Definition runge_kutta_scheme.f90:33

scratch_registry
Defines a registry for storing and requesting temporary objects This can be used when you have a func...
Definition scratch_registry.f90:37

scratch_registry::neko_scratch_registry
type(scratch_registry_t), target, public neko_scratch_registry
Global scratch registry.
Definition scratch_registry.f90:117

ax_product::ax_t
Base type for a matrix-vector product providing .
Definition ax.f90:43

coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:58

euler_res_cpu::euler_res_cpu_t
Definition euler_res_cpu.f90:51

euler_residual::euler_rhs_t
Abstract type to compute rhs.
Definition euler_res.f90:44

field::field_t
Definition field.f90:47

field_list::field_list_t
field_list_t, To be able to group fields together
Definition field_list.f90:15

gather_scatter::gs_t
Definition gather_scatter.f90:68

runge_kutta_time_scheme::runge_kutta_time_scheme_t
Definition runge_kutta_scheme.f90:44