fluid_volflow.f90

Go to the documentation of this file.

! Copyright (c) 2008-2020, UCHICAGO ARGONNE, LLC.
! Copyright (c) 2025-2026, The Neko Authors
!
! The UChicago Argonne, LLC as Operator of Argonne National
! Laboratory holds copyright in the Software. The copyright holder
! reserves all rights except those expressly granted to licensees,
! and U.S. Government license rights.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the disclaimer below.
!
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions and the disclaimer (as noted below)
! in the documentation and/or other materials provided with the
! distribution.
!
! 3. Neither the name of ANL 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
! UCHICAGO ARGONNE, LLC, THE U.S. DEPARTMENT OF
! ENERGY 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.
!
! Additional BSD Notice
! ---------------------
! 1. This notice is required to be provided under our contract with
! the U.S. Department of Energy (DOE). This work was produced at
! Argonne National Laboratory under Contract
! No. DE-AC02-06CH11357 with the DOE.
!
! 2. Neither the United States Government nor UCHICAGO ARGONNE,
! LLC nor any of their employees, makes any warranty,
! express or implied, or assumes any liability or responsibility for the
! accuracy, completeness, or usefulness of any information, apparatus,
! product, or process disclosed, or represents that its use would not
! infringe privately-owned rights.
!
! 3. Also, reference herein to any specific commercial products, process,
! or services by trade name, trademark, manufacturer or otherwise does
! not necessarily constitute or imply its endorsement, recommendation,
! or favoring by the United States Government or UCHICAGO ARGONNE LLC.
! The views and opinions of authors expressed
! herein do not necessarily state or reflect those of the United States
! Government or UCHICAGO ARGONNE, LLC, and shall
! not be used for advertising or product endorsement purposes.
!
module fluid_volflow
  use operators, only : opgrad, cdtp, rotate_cyc
  use num_types, only : rp
  use mathops, only : opchsign
  use krylov, only : ksp_t, ksp_monitor_t
  use precon, only : pc_t
  use dofmap, only : dofmap_t
  use field, only : field_t
  use coefs, only : coef_t
  use time_state, only : time_state_t
  use time_scheme_controller, only : time_scheme_controller_t
  use math, only : copy, glsc2, glmin, glmax, add2, abscmp
  use neko_config, only : neko_bcknd_device
  use device_math, only : device_cfill, device_rzero, device_copy, &
       device_add2, device_add2s2, device_glsc2
  use device_mathops, only : device_opchsign
  use gather_scatter, only : gs_t, gs_op_add
  use json_module, only : json_file
  use json_utils, only : json_get, json_get_or_default, json_get_or_lookup
  use scratch_registry, only : neko_scratch_registry
  use scalar_bc_projector, only : scalar_bc_projector_t
  use vector_bc_projector, only : vector_bc_projector_t
  use ax_product, only : ax_t
  use comm, only : neko_comm, mpi_real_precision, pe_rank
  use mpi_f08, only : mpi_allreduce, mpi_in_place, mpi_sum
  use logger, only : log_size, neko_log
  implicit none
  private
 
  type, public :: fluid_volflow_t
     integer :: flow_dir
     logical :: avflow
     logical :: log = .true.
     real(kind=rp) :: flow_rate
     real(kind=rp) :: dtlag = 0d0
     real(kind=rp) :: bdlag = 0d0 
     type(field_t) :: u_vol, v_vol, w_vol, p_vol
     real(kind=rp) :: domain_length, base_flow
   contains
     procedure, pass(this) :: init => fluid_vol_flow_init
     procedure, pass(this) :: free => fluid_vol_flow_free
     procedure, pass(this) :: adjust => fluid_vol_flow
     procedure, private, pass(this) :: compute => fluid_vol_flow_compute
  end type fluid_volflow_t
 
contains
 
  subroutine fluid_vol_flow_init(this, dm_Xh, params)
    class(fluid_volflow_t), intent(inout) :: this
    type(dofmap_t), target, intent(in) :: dm_Xh
    type(json_file), intent(inout) :: params
    logical average, log_output
    integer :: direction
    real(kind=rp) :: rate
 
    call this%free()
 
    !Initialize vol_flow (if there is a forced volume flow)
    call json_get_or_lookup(params, 'case.fluid.flow_rate_force.direction', &
         direction)
    call json_get_or_lookup(params, 'case.fluid.flow_rate_force.value', rate)
    call json_get(params, 'case.fluid.flow_rate_force.use_averaged_flow',&
         average)
    call json_get_or_default(params, 'case.fluid.flow_rate_force.log', &
         log_output, .true.)
 
    this%flow_dir = direction
    this%avflow = average
    this%log = log_output
    this%flow_rate = rate
 
    if (this%flow_dir .ne. 0) then
       call this%u_vol%init(dm_xh, 'u_vol')
       call this%v_vol%init(dm_xh, 'v_vol')
       call this%w_vol%init(dm_xh, 'w_vol')
       call this%p_vol%init(dm_xh, 'p_vol')
    end if
 
  end subroutine fluid_vol_flow_init
 
  subroutine fluid_vol_flow_free(this)
    class(fluid_volflow_t), intent(inout) :: this
 
    call this%u_vol%free()
    call this%v_vol%free()
    call this%w_vol%free()
    call this%p_vol%free()
 
  end subroutine fluid_vol_flow_free
 
  subroutine fluid_vol_flow_compute(this, u_res, v_res, w_res, p_res, &
       ext_bdf, gs_Xh, c_Xh, rho, mu, bd, dt, &
       dp_projector, vel_projector, &
       Ax_vel, Ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, &
       vel_max_iter)
    class(fluid_volflow_t), intent(inout) :: this
    type(field_t), intent(inout) :: u_res, v_res, w_res, p_res
    type(coef_t), intent(inout) :: c_Xh
    type(gs_t), intent(inout) :: gs_Xh
    type(time_scheme_controller_t), intent(in) :: ext_bdf
    type(scalar_bc_projector_t), intent(inout) :: dp_projector
    class(vector_bc_projector_t), intent(inout) :: vel_projector
    class(ax_t), intent(in) :: Ax_vel
    class(ax_t), intent(in) :: Ax_prs
    class(ksp_t), intent(inout) :: ksp_prs, ksp_vel
    class(pc_t), intent(inout) :: pc_prs, pc_vel
    real(kind=rp), intent(in) :: bd
    real(kind=rp), intent(in) :: rho, dt
    type(field_t) :: mu
    integer, intent(in) :: vel_max_iter, prs_max_iter
    integer :: n, i
    real(kind=rp) :: xlmin, xlmax
    real(kind=rp) :: ylmin, ylmax
    real(kind=rp) :: zlmin, zlmax
    type(ksp_monitor_t) :: ksp_results(4)
    type(field_t), pointer :: ta1, ta2, ta3
    integer :: temp_indices(3)
 
    call neko_scratch_registry%request_field(ta1, temp_indices(1), .false.)
    call neko_scratch_registry%request_field(ta2, temp_indices(2), .false.)
    call neko_scratch_registry%request_field(ta3, temp_indices(3), .false.)
 
 
    associate(msh => c_xh%msh, p_vol => this%p_vol, &
         u_vol => this%u_vol, v_vol => this%v_vol, w_vol => this%w_vol)
 
      n = c_xh%dof%size()
      xlmin = glmin(c_xh%dof%x, n)
      xlmax = glmax(c_xh%dof%x, n)
      ylmin = glmin(c_xh%dof%y, n) !  for Y!
      ylmax = glmax(c_xh%dof%y, n)
      zlmin = glmin(c_xh%dof%z, n) !  for Z!
      zlmax = glmax(c_xh%dof%z, n)
      if (this%flow_dir .eq. 1) then
         this%domain_length = xlmax - xlmin
      end if
      if (this%flow_dir .eq. 2) then
         this%domain_length = ylmax - ylmin
      end if
      if (this%flow_dir .eq. 3) then
         this%domain_length = zlmax - zlmin
      end if
 
      if (neko_bcknd_device .eq. 1) then
         call device_cfill(c_xh%h1_d, 1.0_rp/rho, n)
         call device_rzero(c_xh%h2_d, n)
      else
         !OCL NORECURRENCE, NOVREC, NOALIAS
         !DIR$ CONCURRENT
         !DIR$ IVDEP
         !GCC$ ivdep
         !$omp parallel do
         do i = 1, n
            c_xh%h1(i,1,1,1) = 1.0_rp / rho
            c_xh%h2(i,1,1,1) = 0.0_rp
         end do
         !$omp end parallel do
      end if
      c_xh%ifh2 = .false.
 
      !   Compute pressure
 
      if (this%flow_dir .eq. 1) then
         call cdtp(p_res%x, c_xh%h1, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
      end if
 
      if (this%flow_dir .eq. 2) then
         call cdtp(p_res%x, c_xh%h1, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
      end if
 
      if (this%flow_dir .eq. 3) then
         call cdtp(p_res%x, c_xh%h1, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
      end if
 
      call gs_xh%op(p_res, gs_op_add)
      call dp_projector%apply(p_res%x, n)
      call pc_prs%update()
      ksp_results(1) = ksp_prs%solve(ax_prs, p_vol, p_res%x, n, &
           c_xh, dp_projector, gs_xh, prs_max_iter)
 
      !   Compute velocity
 
      call opgrad(u_res%x, v_res%x, w_res%x, p_vol%x, c_xh)
 
      if (neko_bcknd_device .eq. 1) then
         call device_opchsign(u_res%x_d, v_res%x_d, w_res%x_d, msh%gdim, n)
         call device_copy(ta1%x_d, c_xh%B_d, n)
         call device_copy(ta2%x_d, c_xh%B_d, n)
         call device_copy(ta3%x_d, c_xh%B_d, n)
      else
         call opchsign(u_res%x, v_res%x, w_res%x, msh%gdim, n)
         call copy(ta1%x, c_xh%B, n)
         call copy(ta2%x, c_xh%B, n)
         call copy(ta3%x, c_xh%B, n)
      end if
      call vel_projector%apply(ta1%x, ta2%x, ta3%x, n)
 
      ! add forcing
 
      if (neko_bcknd_device .eq. 1) then
         if (this%flow_dir .eq. 1) then
            call device_add2(u_res%x_d, ta1%x_d, n)
         else if (this%flow_dir .eq. 2) then
            call device_add2(v_res%x_d, ta2%x_d, n)
         else if (this%flow_dir .eq. 3) then
            call device_add2(w_res%x_d, ta3%x_d, n)
         end if
      else
         if (this%flow_dir .eq. 1) then
            call add2(u_res%x, ta1%x, n)
         else if (this%flow_dir .eq. 2) then
            call add2(v_res%x, ta2%x, n)
         else if (this%flow_dir .eq. 3) then
            call add2(w_res%x, ta3%x, n)
         end if
      end if
 
      if (neko_bcknd_device .eq. 1) then
         call device_copy(c_xh%h1_d, mu%x_d, n)
         call device_cfill(c_xh%h2_d, rho * (bd / dt), n)
      else
         call copy(c_xh%h1, mu%x, n)
         c_xh%h2 = rho * (bd / dt)
      end if
      c_xh%ifh2 = .true.
 
      call rotate_cyc(u_res, v_res, w_res, 1, c_xh)
      call gs_xh%op(u_res%x, v_res%x, w_res%x, n, gs_op_add)
      call rotate_cyc(u_res, v_res, w_res, 0, c_xh)
 
      call vel_projector%apply(u_res%x, v_res%x, w_res%x, n)
      call pc_vel%update()
 
      ksp_results(2:4) = ksp_vel%solve_coupled(ax_vel, &
           u_vol, v_vol, w_vol, &
           u_res%x, v_res%x, w_res%x, &
           n, c_xh, &
           vel_projector, &
           gs_xh, vel_max_iter)
 
      if (neko_bcknd_device .eq. 1) then
         if (this%flow_dir .eq. 1) then
            this%base_flow = &
                 device_glsc2(u_vol%x_d, c_xh%B_d, n) / this%domain_length
         end if
 
         if (this%flow_dir .eq. 2) then
            this%base_flow = &
                 device_glsc2(v_vol%x_d, c_xh%B_d, n) / this%domain_length
         end if
 
         if (this%flow_dir .eq. 3) then
            this%base_flow = &
                 device_glsc2(w_vol%x_d, c_xh%B_d, n) / this%domain_length
         end if
      else
         if (this%flow_dir .eq. 1) then
            this%base_flow = glsc2(u_vol%x, c_xh%B, n) / this%domain_length
         end if
 
         if (this%flow_dir .eq. 2) then
            this%base_flow = glsc2(v_vol%x, c_xh%B, n) / this%domain_length
         end if
 
         if (this%flow_dir .eq. 3) then
            this%base_flow = glsc2(w_vol%x, c_xh%B, n) / this%domain_length
         end if
      end if
    end associate
 
    call neko_scratch_registry%relinquish_field(temp_indices)
  end subroutine fluid_vol_flow_compute
 
  subroutine fluid_vol_flow(this, u, v, w, p, u_res, v_res, w_res, p_res, &
       c_Xh, gs_Xh, ext_bdf, rho, mu, dt, time, &
       dp_projector, vel_projector, &
       Ax_vel, Ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, &
       vel_max_iter)
 
    class(fluid_volflow_t), intent(inout) :: this
    type(field_t), intent(inout) :: u, v, w, p
    type(field_t), intent(inout) :: u_res, v_res, w_res, p_res
    type(coef_t), intent(inout) :: c_Xh
    type(gs_t), intent(inout) :: gs_Xh
    type(time_scheme_controller_t), intent(in) :: ext_bdf
    type(time_state_t), intent(in) :: time
    real(kind=rp), intent(in) :: rho, dt
    type(field_t) :: mu
    type(scalar_bc_projector_t), intent(inout) :: dp_projector
    class(vector_bc_projector_t), intent(inout) :: vel_projector
    class(ax_t), intent(in) :: Ax_vel
    class(ax_t), intent(in) :: Ax_prs
    class(ksp_t), intent(inout) :: ksp_prs, ksp_vel
    class(pc_t), intent(inout) :: pc_prs, pc_vel
    integer, intent(in) :: prs_max_iter, vel_max_iter
    real(kind=rp) :: ifcomp, flow_rate, xsec
    real(kind=rp) :: current_flow, delta_flow, scale
    integer :: n, ierr, i
    character(len=5) :: flow_dir_label
    character(len=12) :: step_str
 
    character(len=200) :: log_buf
 
    associate(u_vol => this%u_vol, v_vol => this%v_vol, &
         w_vol => this%w_vol, p_vol => this%p_vol)
 
      n = c_xh%dof%size()
 
      ! If either dt or the backwards difference coefficient change,
      ! then recompute base flow solution corresponding to unit forcing:
 
      ifcomp = 0.0_rp
 
      if ((.not. abscmp(dt, this%dtlag)) .or. &
           (.not. abscmp(ext_bdf%diffusion_coeffs%x(1), this%bdlag))) then
         ifcomp = 1.0_rp
      end if
 
      this%dtlag = dt
      this%bdlag = ext_bdf%diffusion_coeffs%x(1)
 
      call mpi_allreduce(mpi_in_place, ifcomp, 1, &
           mpi_real_precision, mpi_sum, neko_comm, ierr)
 
      if (ifcomp .gt. 0d0) then
         call this%compute(u_res, v_res, w_res, p_res, &
              ext_bdf, gs_xh, c_xh, rho, mu, ext_bdf%diffusion_coeffs%x(1), dt, &
              dp_projector, vel_projector, &
              ax_vel, ax_prs, ksp_prs, ksp_vel, pc_prs, pc_vel, prs_max_iter, &
              vel_max_iter)
      end if
 
      if (neko_bcknd_device .eq. 1) then
         if (this%flow_dir .eq. 1) then
            current_flow = &
                 device_glsc2(u%x_d, c_xh%B_d, n) / this%domain_length ! for X
         else if (this%flow_dir .eq. 2) then
            current_flow = &
                 device_glsc2(v%x_d, c_xh%B_d, n) / this%domain_length ! for Y
         else if (this%flow_dir .eq. 3) then
            current_flow = &
                 device_glsc2(w%x_d, c_xh%B_d, n) / this%domain_length ! for Z
         end if
      else
         if (this%flow_dir .eq. 1) then
            current_flow = glsc2(u%x, c_xh%B, n) / this%domain_length ! for X
         else if (this%flow_dir .eq. 2) then
            current_flow = glsc2(v%x, c_xh%B, n) / this%domain_length ! for Y
         else if (this%flow_dir .eq. 3) then
            current_flow = glsc2(w%x, c_xh%B, n) / this%domain_length ! for Z
         end if
      end if
 
      if (this%avflow) then
         xsec = c_xh%volume / this%domain_length
         flow_rate = this%flow_rate*xsec
      else
         flow_rate = this%flow_rate
      end if
 
      delta_flow = flow_rate - current_flow
      scale = delta_flow / this%base_flow
 
      if (this%log .and. pe_rank .eq. 0) then
         if (this%flow_dir .eq. 1) then
            flow_dir_label = '    x'
         else if (this%flow_dir .eq. 2) then
            flow_dir_label = '    y'
         else if (this%flow_dir .eq. 3) then
            flow_dir_label = '    z'
         end if
         write(step_str, '(I12)') time%tstep
         step_str = adjustl(step_str)
         write(log_buf, '(A,A3,A5,1X,5A18)') &
              'Flow rate   ', ' | ', 'Dir.:', 'Time:', 'Scale:', 'Rate:', &
              'Current:', 'Base:'
         call neko_log%message(log_buf)
         write(log_buf, '(A12,A3,A5,1X,5E18.9)') step_str, ' | ', &
              flow_dir_label, time%t, scale, flow_rate, current_flow, &
              this%base_flow
         call neko_log%message(log_buf)
      end if
 
      if (neko_bcknd_device .eq. 1) then
         call device_add2s2(u%x_d, u_vol%x_d, scale, n)
         call device_add2s2(v%x_d, v_vol%x_d, scale, n)
         call device_add2s2(w%x_d, w_vol%x_d, scale, n)
         call device_add2s2(p%x_d, p_vol%x_d, scale, n)
      else
         !OCL NORECURRENCE, NOVREC, NOALIAS
         !DIR$ CONCURRENT
         !DIR$ IVDEP
         !GCC$ ivdep
         !$omp parallel do
         do i = 1, n
            u%x(i,1,1,1) = u%x(i,1,1,1) + scale * u_vol%x(i,1,1,1)
            v%x(i,1,1,1) = v%x(i,1,1,1) + scale * v_vol%x(i,1,1,1)
            w%x(i,1,1,1) = w%x(i,1,1,1) + scale * w_vol%x(i,1,1,1)
            p%x(i,1,1,1) = p%x(i,1,1,1) + scale * p_vol%x(i,1,1,1)
         end do
         !$omp end parallel do
      end if
    end associate
 
  end subroutine fluid_vol_flow
 
end module fluid_volflow