ale_routines_cpu.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.
!
! Routines for expensive ALE calculations on CPU
module ale_routines_cpu
  use num_types, only : rp
  use field, only : field_t
  use coefs, only : coef_t
  use math, only : cfill, glimax, rzero
  use field_series, only : field_series_t
  use time_state, only : time_state_t
  use ab_time_scheme, only : ab_time_scheme_t
  use utils, only : neko_error
  use comm, only : neko_comm
  use zero_dirichlet, only : zero_dirichlet_t
  use mesh, only : mesh_t
  use gather_scatter, only : gs_op_min
  use mpi_f08, only : mpi_wtime, mpi_barrier
  use logger, only : neko_log, log_size
  use ale_rigid_kinematics, only : ale_config_t, body_kinematics_t
  implicit none
  private
 
  public :: add_kinematics_to_mesh_velocity_cpu
  public :: update_ale_mesh_cpu
  public :: compute_cheap_dist_v2_cpu
 
contains
 
  subroutine compute_cheap_dist_cpu(d, coef, msh, zone_indices)
    real(kind=rp), intent(inout), target :: d(:)
    type(coef_t), intent(in) :: coef
    type(mesh_t), intent(in) :: msh
    type(zero_dirichlet_t) :: bc_wall
    integer, intent(in) :: zone_indices(:)
    real(kind=rp), pointer :: d4(:, :, :, :)
    integer :: i, j, k, e, n
    integer :: ipass, nchange, max_pass
    integer :: ii, jj, kk, i0, i1, j0, j1, k0, k1
    integer :: lx, ly, lz, nel, z_idx
    integer :: m, idx
    real(kind=rp) :: dtmp, x1, y1, z1, x2, y2, z2
    integer :: change_vec(1)
    logical :: done
 
    lx = coef%dof%Xh%lx
    ly = coef%dof%Xh%ly
    lz = coef%dof%Xh%lz
    nel = msh%nelv
    n = size(d)
    d4(1:lx, 1:ly, 1:lz, 1:nel) => d
    max_pass = 10000
 
    !d = huge(0.0_rp)
    call cfill(d, huge(0.0_rp), n)
 
    if (size(zone_indices) .gt. 0) then
       call bc_wall%init_from_components(coef)
       do k = 1, size(zone_indices)
          z_idx = zone_indices(k)
          call bc_wall%mark_zone(msh%labeled_zones(z_idx))
       end do
       call bc_wall%finalize()
       m = bc_wall%msk(0)
       do i = 1, m
          idx = bc_wall%msk(i)
          d(idx) = 0.0_rp
       end do
       call bc_wall%free()
    end if
 
    ipass = 1
    done = .false.
    do while ((ipass .le. max_pass) .and. .not. done)
       nchange = 0
       do e = 1, nel
          do k = 1, lz
             do j = 1, ly
                do i = 1, lx
                   x1 = coef%dof%x(i, j, k, e)
                   y1 = coef%dof%y(i, j, k, e)
                   z1 = coef%dof%z(i, j, k, e)
                   i0 = max(1, i - 1)
                   i1 = min(lx, i + 1)
                   j0 = max(1, j - 1)
                   j1 = min(ly, j + 1)
                   k0 = max(1, k - 1)
                   k1 = min(lz, k + 1)
                   do kk = k0, k1
                      do jj = j0, j1
                         do ii = i0, i1
                            if ((ii .eq. i) .and. (jj .eq. j) .and. &
                                 (kk .eq. k)) cycle
                            x2 = coef%dof%x(ii, jj, kk, e)
                            y2 = coef%dof%y(ii, jj, kk, e)
                            z2 = coef%dof%z(ii, jj, kk, e)
                            dtmp = d4(ii, jj, kk, e) + &
                                 sqrt((x1 - x2)**2 + &
                                 (y1 - y2)**2 + (z1 - z2)**2)
                            if (dtmp .lt. d4(i, j, k, e)) then
                               d4(i, j, k, e) = dtmp
                               nchange = nchange + 1
                            end if
                         end do
                      end do
                   end do
                end do
             end do
          end do
       end do
       call coef%gs_h%gs_op_vector(d, n, gs_op_min)
       change_vec(1) = nchange
       if (glimax(change_vec, 1) .eq. 0) done = .true.
       ipass = ipass + 1
    end do
  end subroutine compute_cheap_dist_cpu
 
  subroutine compute_cheap_dist_v2_cpu(dist_field, coef, msh, zone_indices)
    type(field_t), intent(inout) :: dist_field
    type(coef_t), intent(in) :: coef
    type(mesh_t), intent(in) :: msh
    integer, intent(in) :: zone_indices(:)
    type(zero_dirichlet_t) :: bc_wall
    integer :: i, j, k, e, n
    integer :: ipass, nchange, max_pass
    integer :: ii, jj, kk, i0, i1, j0, j1, k0, k1
    integer :: lx, ly, lz, nel, z_idx
    integer :: m, idx, iter, local_iters
    real(kind=rp) :: dtmp, x1, y1, z1, x2, y2, z2
    integer :: change_vec(1)
    logical :: done, changed_local, element_changed_ever
    character(len=LOG_SIZE) :: log_buf
 
    lx = coef%dof%Xh%lx
    ly = coef%dof%Xh%ly
    lz = coef%dof%Xh%lz
    nel = msh%nelv
    n = coef%dof%size()
    max_pass = 10000
 
    ! Limit for worst case scenario such that all nodes can propagate
    ! their values across the element before triggering an MPI call.
    local_iters = lx + ly + lz
 
    call cfill(dist_field%x, huge(0.0_rp), n)
 
    if (size(zone_indices) .gt. 0) then
       call bc_wall%init_from_components(coef)
       do k = 1, size(zone_indices)
          z_idx = zone_indices(k)
          call bc_wall%mark_zone(msh%labeled_zones(z_idx))
       end do
       call bc_wall%finalize()
       m = bc_wall%msk(0)
       do i = 1, m
          idx = bc_wall%msk(i)
          dist_field%x(idx, 1, 1, 1) = 0.0_rp
       end do
       call bc_wall%free()
    end if
 
    ipass = 1
    done = .false.
    do while ((ipass .le. max_pass) .and. .not. done)
       nchange = 0
 
       do e = 1, nel
          iter = 1
          element_changed_ever = .false.
          changed_local = .true.
          do while (changed_local .and. (iter .le. local_iters))
 
             changed_local = .false.
             do k = 1, lz
                do j = 1, ly
                   do i = 1, lx
                      x1 = coef%dof%x(i, j, k, e)
                      y1 = coef%dof%y(i, j, k, e)
                      z1 = coef%dof%z(i, j, k, e)
                      i0 = max(1, i - 1)
                      i1 = min(lx, i + 1)
                      j0 = max(1, j - 1)
                      j1 = min(ly, j + 1)
                      k0 = max(1, k - 1)
                      k1 = min(lz, k + 1)
                      do kk = k0, k1
                         do jj = j0, j1
                            do ii = i0, i1
                               if ((ii .eq. i) .and. (jj .eq. j) .and. &
                                    (kk .eq. k)) cycle
 
                               x2 = coef%dof%x(ii, jj, kk, e)
                               y2 = coef%dof%y(ii, jj, kk, e)
                               z2 = coef%dof%z(ii, jj, kk, e)
 
                               dtmp = dist_field%x(ii, jj, kk, e) + &
                                    sqrt((x1 - x2)**2 + &
                                    (y1 - y2)**2 + (z1 - z2)**2)
 
                               if (dtmp .lt. dist_field%x(i, j, k, e)) then
                                  dist_field%x(i, j, k, e) = dtmp
                                  changed_local = .true.
                               end if
 
                            end do
                         end do
                      end do
                   end do
                end do
             end do
             if (changed_local) element_changed_ever = .true.
             iter = iter + 1
          end do
 
          if (element_changed_ever) nchange = nchange + 1
       end do
 
       call coef%gs_h%gs_op_vector(dist_field%x, n, gs_op_min)
       change_vec(1) = nchange
 
       if (glimax(change_vec, 1) .eq. 0) done = .true.
       ipass = ipass + 1
    end do
 
    if (done) then
       write(log_buf, '(A, I0, A)') "    converged in: ", ipass, " passes"
       call neko_log%message(log_buf)
    else
       write(log_buf, '(A, I0, A)') "    reached max passes: ", ipass, &
            " without convergence"
       call neko_log%message(log_buf)
    end if
  end subroutine compute_cheap_dist_v2_cpu
 
  subroutine add_kinematics_to_mesh_velocity_cpu(wx, wy, wz, &
       x_ref, y_ref, z_ref, phi, coef, kinematics, rot_mat, inital_pivot_loc)
    type(field_t), intent(inout) :: wx, wy, wz
    type(field_t), intent(in) :: x_ref, y_ref, z_ref
    type(field_t), intent(in) :: phi
    type(coef_t), intent(in) :: coef
    type(body_kinematics_t), intent(in) :: kinematics
    real(kind=rp), intent(in) :: inital_pivot_loc(3)
    real(kind=rp), intent(in) :: rot_mat(3,3)
    integer :: i, n
    real(kind=rp) :: rx, ry, rz
    real(kind=rp) :: v_tan_x, v_tan_y, v_tan_z
    real(kind=rp) :: dx_ref, dy_ref, dz_ref
    real(kind=rp) :: rx_target, ry_target, rz_target
    real(kind=rp) :: p_val
    n = phi%dof%size()
    associate(x => coef%dof%x, y => coef%dof%y, z => coef%dof%z)
      do concurrent(i = 1:n)
         ! for points on the wall (phi=1) we do this to avoid any
         ! numeric error due to computation of rotation matrix.
         ! It ensures, the walls are always where they need to be!
         if ( abs(phi%x(i, 1, 1, 1) - 1.0_rp) .lt. 1e-6_rp ) then
 
            rx = x(i, 1, 1, 1) - kinematics%center(1)
            ry = y(i, 1, 1, 1) - kinematics%center(2)
            rz = z(i, 1, 1, 1) - kinematics%center(3)
            v_tan_x = kinematics%vel_ang(2) * rz - kinematics%vel_ang(3) * ry
            v_tan_y = kinematics%vel_ang(3) * rx - kinematics%vel_ang(1) * rz
            v_tan_z = kinematics%vel_ang(1) * ry - kinematics%vel_ang(2) * rx
            wx%x(i, 1, 1, 1) = wx%x(i, 1, 1, 1) + &
                  (kinematics%vel_trans(1) + v_tan_x) * phi%x(i, 1, 1, 1)
            wy%x(i, 1, 1, 1) = wy%x(i, 1, 1, 1) + &
                  (kinematics%vel_trans(2) + v_tan_y) * phi%x(i, 1, 1, 1)
            wz%x(i, 1, 1, 1) = wz%x(i, 1, 1, 1) + &
                  (kinematics%vel_trans(3) + v_tan_z) * phi%x(i, 1, 1, 1)
         else
            ! For other points we do this to avoid the time-dependnent
            ! drift in some special cases, which happens due to the nature of
            ! the ODE that we integrate above!
            dx_ref = x_ref%x(i, 1, 1, 1) - inital_pivot_loc(1)
            dy_ref = y_ref%x(i, 1, 1, 1) - inital_pivot_loc(2)
            dz_ref = z_ref%x(i, 1, 1, 1) - inital_pivot_loc(3)
 
            ! Rotate to find the "ghost" target vector
            ! Apply the rotation matrix R(t) to the reference vector
            rx_target = rot_mat(1,1)*dx_ref + rot_mat(1,2)*dy_ref + &
                 rot_mat(1,3)*dz_ref
            ry_target = rot_mat(2,1)*dx_ref + rot_mat(2,2)*dy_ref + &
                 rot_mat(2,3)*dz_ref
            rz_target = rot_mat(3,1)*dx_ref + rot_mat(3,2)*dy_ref + &
                 rot_mat(3,3)*dz_ref
 
            ! Calculate tangential velocity at the ghost target
            ! v_tan = Omega \corss R_target
            v_tan_x = kinematics%vel_ang(2) * rz_target - &
                 kinematics%vel_ang(3) * ry_target
            v_tan_y = kinematics%vel_ang(3) * rx_target - &
                 kinematics%vel_ang(1) * rz_target
            v_tan_z = kinematics%vel_ang(1) * ry_target - &
                 kinematics%vel_ang(2) * rx_target
 
            p_val = phi%x(i, 1, 1, 1)
 
            ! Total Mesh Velocity
            wx%x(i, 1, 1, 1) = wx%x(i, 1, 1, 1) + &
                 (kinematics%vel_trans(1) + v_tan_x) * p_val
            wy%x(i, 1, 1, 1) = wy%x(i, 1, 1, 1) + &
                 (kinematics%vel_trans(2) + v_tan_y) * p_val
            wz%x(i, 1, 1, 1) = wz%x(i, 1, 1, 1) + &
                 (kinematics%vel_trans(3) + v_tan_z) * p_val
 
         end if
      end do
    end associate
  end subroutine add_kinematics_to_mesh_velocity_cpu
 
  subroutine update_ale_mesh_cpu(c_Xh, wm_x, wm_y, wm_z, wm_x_lag, wm_y_lag, &
       wm_z_lag, time, nadv, scheme_type)
    type(coef_t), intent(inout) :: c_xh
    type(field_t), intent(in) :: wm_x, wm_y, wm_z
    type(field_series_t), intent(in) :: wm_x_lag, wm_y_lag, wm_z_lag
    type(time_state_t), intent(in) :: time
    type(ab_time_scheme_t) :: ab_scheme_obj
    integer, intent(in) :: nadv
    integer :: i, j, n
    character(len=*), intent(in) :: scheme_type
    real(kind=rp) :: ab_coeffs(4), dt_history(10)
 
    call rzero(ab_coeffs, 4)
    if (trim(scheme_type) .eq. 'ab') then
       dt_history(1) = time%dt
       dt_history(2) = time%dtlag(1)
       dt_history(3) = time%dtlag(2)
       call ab_scheme_obj%compute_coeffs(ab_coeffs, dt_history, nadv)
    else
       call neko_error("ALE: Unknown mesh time-integration scheme")
    end if
 
    n = c_xh%dof%size()
 
    do concurrent(i = 1:n)
       c_xh%dof%x(i, 1, 1, 1) = c_xh%dof%x(i, 1, 1, 1) + &
            time%dt * ab_coeffs(1) * wm_x%x(i, 1, 1, 1)
       c_xh%dof%y(i, 1, 1, 1) = c_xh%dof%y(i, 1, 1, 1) + &
            time%dt * ab_coeffs(1) * wm_y%x(i, 1, 1, 1)
       c_xh%dof%z(i, 1, 1, 1) = c_xh%dof%z(i, 1, 1, 1) + &
            time%dt * ab_coeffs(1) * wm_z%x(i, 1, 1, 1)
 
       do j = 2, nadv
          c_xh%dof%x(i, 1, 1, 1) = c_xh%dof%x(i, 1, 1, 1) + &
               time%dt * ab_coeffs(j) * wm_x_lag%lf(j - 1)%x(i, 1, 1, 1)
          c_xh%dof%y(i, 1, 1, 1) = c_xh%dof%y(i, 1, 1, 1) + &
               time%dt * ab_coeffs(j) * wm_y_lag%lf(j - 1)%x(i, 1, 1, 1)
          c_xh%dof%z(i, 1, 1, 1) = c_xh%dof%z(i, 1, 1, 1) + &
               time%dt * ab_coeffs(j) * wm_z_lag%lf(j - 1)%x(i, 1, 1, 1)
       end do
 
    end do
  end subroutine update_ale_mesh_cpu
 
end module ale_routines_cpu