fld_file_data.f90

Go to the documentation of this file.

 
module fld_file_data
  use num_types, only : rp, dp
  use math, only : cmult, add2
  use vector, only : vector_t, vector_ptr_t
  use interpolation, only : interpolator_t
  use field, only : field_t
  use field_list, only : field_list_t
  use logger, only : neko_log, log_size
  use device, only : host_to_device
  use dofmap, only : dofmap_t
  use json_utils, only : json_get_or_default
  use json_module, only : json_file
  use space, only : space_t, gll
  use global_interpolation, only : global_interpolation_t, &
       glob_interp_pad, glob_interp_tol
  use utils, only : neko_error, neko_fname_len, extract_fld_file_index
  use mesh, only : mesh_t
  implicit none
  private
 
  type, public :: fld_file_data_t
     type(vector_t) :: x
     type(vector_t) :: y
     type(vector_t) :: z
     type(vector_t) :: u
     type(vector_t) :: v
     type(vector_t) :: w
     type(vector_t) :: p
     type(vector_t) :: t
     integer, allocatable :: idx(:)
     type(vector_t), allocatable :: s(:)
     integer :: gdim
     integer :: n_scalars = 0
     real(kind=rp) :: time = 0.0 
     integer :: glb_nelv = 0
     integer :: nelv = 0
     integer :: offset_el = 0
     integer :: lx = 0
     integer :: ly = 0
     integer :: lz = 0
     integer :: t_counter = 0
     ! meta file information (if any)
     integer :: meta_nsamples = 0
     integer :: meta_start_counter = 0
     character(len=1024) :: fld_series_fname
   contains
     procedure, pass(this) :: init => fld_file_data_init
     procedure, pass(this) :: free => fld_file_data_free
     procedure, pass(this) :: scale => fld_file_data_scale
     procedure, pass(this) :: add => fld_file_data_add
     procedure, pass(this) :: size => fld_file_data_size
     procedure, pass(this) :: get_list => fld_file_data_get_list
     procedure, pass(this) :: init_same => fld_file_data_init_same
     procedure, pass(this) :: init_n_fields => fld_file_data_init_n_fields
     procedure, pass(this) :: generate_interpolator => &
          fld_file_data_generate_interpolator
     procedure, pass(this) :: import_fields => fld_file_data_import_fields
  end type fld_file_data_t
 
contains
 
  subroutine fld_file_data_import_fields(this, u, v, w, p, t, &
       s_target_list, s_index_list, interpolate, tolerance, padding, &
       global_interp_subdict)
    class(fld_file_data_t), intent(inout) :: this
    type(field_t), pointer, intent(inout), optional :: u,v,w,p,t
    type(field_list_t), intent(inout), optional :: s_target_list
    integer, intent(in), optional :: s_index_list(:)
    logical, intent(in), optional :: interpolate
    real(kind=dp), intent(in), optional :: tolerance
    real(kind=dp), intent(in), optional :: padding
    type(json_file), intent(inout), optional :: global_interp_subdict
 
    integer :: i
 
    ! ---- For the mesh to mesh interpolation
    logical :: mesh_mismatch
    logical :: interpolate_
    type(global_interpolation_t) :: global_interp
    type(dofmap_t), pointer :: dof
    type(mesh_t) , pointer :: msh
    ! -----
 
    ! ---- For space to space interpolation
    type(space_t) :: prev_Xh
    type(space_t) , pointer :: Xh
    type(interpolator_t) :: space_interp
    ! ----
 
    character(len=LOG_SIZE) :: log_buf
 
    ! ---- Default values
    interpolate_ = .false.
    if (present(interpolate)) interpolate_ = interpolate
    ! ----
 
    !
    ! Handle the passing of arguments and pointers
    !
    dof => null()
    msh => null()
    xh => null()
 
    if (present(u)) then
       dof => u%dof; msh => u%msh; xh => u%Xh
    else if (present(v)) then
       dof => v%dof; msh => v%msh; xh => v%Xh
    else if (present(w)) then
       dof => w%dof; msh => w%msh; xh => w%Xh
    else if (present(p)) then
       dof => p%dof; msh => p%msh; xh => p%Xh
    else if (present(t)) then
       dof => t%dof; msh => t%msh; xh => t%Xh
    else if (present(s_target_list)) then
       if (s_target_list%size() .eq. 0) then
          call neko_error("Scalar target list is empty")
       else
          dof => s_target_list%items(1)%ptr%dof
          msh => s_target_list%items(1)%ptr%msh
          xh => s_target_list%items(1)%ptr%Xh
       end if
    else
       call neko_error("At least one field must be passed")
    end if
 
    !
    ! Check that the data in the fld file matches the current case.
    ! Note that this is a safeguard and there are corner cases where
    ! two different meshes have the same dimension and same # of elements
    ! but this should be enough to cover obvious cases.
    !
    mesh_mismatch = (this%glb_nelv .ne. msh%glb_nelv .or. &
         this%gdim .ne. msh%gdim)
 
    if (mesh_mismatch .and. .not. interpolate_) then
       call neko_error("The fld file must match the current mesh! &
       &Use 'interpolate': 'true' to enable interpolation.")
    else if (.not. mesh_mismatch .and. interpolate_) then
       call neko_log%warning("You have activated interpolation but you might &
       &still be using the same mesh.")
    end if
 
    !
    ! Handling of interpolation and I/O
    !
    if (interpolate_) then
 
       ! Throw error if dof or msh are not specified
       ! This should never be thrown, but just in case.
       if (.not. associated(dof) .or. .not. associated(msh)) then
          call neko_error("both dof and msh must be associated")
       end if
 
       ! Generates an interpolator object and performs the point search
       ! NOTE: optional arguments are managed in the subroutine, but if both
       ! are provided the JSON subdict takes precedence
       call this%generate_interpolator(global_interp, dof, msh, &
            tolerance = tolerance, padding = padding, &
            global_interp_subdict = global_interp_subdict)
 
       ! Evaluate all the fields
       if (present(u)) call global_interp%evaluate(u%x(:,1,1,1), this%u%x, &
            on_host = .false.)
       if (present(v)) call global_interp%evaluate(v%x(:,1,1,1), this%v%x, &
            on_host = .false.)
       if (present(w)) call global_interp%evaluate(w%x(:,1,1,1), this%w%x, &
            on_host = .false.)
       if (present(p)) call global_interp%evaluate(p%x(:,1,1,1), this%p%x, &
            on_host = .false.)
       if (present(t)) call global_interp%evaluate(t%x(:,1,1,1), this%t%x, &
            on_host = .false.)
       if (present(s_target_list)) then
 
          ! If the index list exists, use it as a "mask"
          if (present(s_index_list)) then
             do i = 1, size(s_index_list)
                ! Take care that if we set i=0 we want temperature
                if (s_index_list(i) .eq. 0) then
                   call global_interp%evaluate(s_target_list%x(i), &
                        this%t%x, on_host = .false.)
                else
                   ! For scalar fields, require indices in 1:this%n_scalars
                   if (s_index_list(i) < 1 .or. &
                        s_index_list(i) > this%n_scalars) then
                      call neko_error("s_index_list entry out of bounds")
                   end if
                   call global_interp%evaluate(s_target_list%x(i), &
                        this%s(s_index_list(i))%x, on_host = .false.)
                end if
             end do
 
             ! otherwise, just copy element-to-element
          else
             do i = 1, s_target_list%size()
                call global_interp%evaluate(s_target_list%x(i), this%s(i)%x, &
                     on_host = .false.)
             end do
          end if ! present s_index_list
       end if ! present s_tgt
 
       call global_interp%free()
 
    else ! No interpolation, but potentially just from different spaces
 
       ! throw an error is the space is not passed
       if (.not. associated(xh)) call neko_error("Xh is not associated")
 
       ! Build a space_t object from the data in the fld file
       call prev_xh%init(gll, this%lx, this%ly, this%lz)
       call space_interp%init(xh, prev_xh)
 
       ! Do the space-to-space interpolation
       if (present(u)) call space_interp%map(u%x, this%u%x, this%nelv, xh)
       if (present(v)) call space_interp%map(v%x, this%v%x, this%nelv, xh)
       if (present(w)) call space_interp%map(w%x, this%w%x, this%nelv, xh)
       if (present(p)) call space_interp%map(p%x, this%p%x, this%nelv, xh)
       if (present(t)) call space_interp%map(t%x, this%t%x, this%nelv, xh)
       if (present(s_target_list)) then
 
          ! If the index list exists, use it as a "mask"
          if (present(s_index_list)) then
             do i = 1, size(s_index_list)
 
                ! 0 means we want temperature
                if (s_index_list(i) .eq. 0) then
                   call space_interp%map(s_target_list%x(i), &
                        this%t%x, this%nelv, xh)
                else
                   call space_interp%map(s_target_list%x(i), &
                        this%s(s_index_list(i))%x, this%nelv, xh)
                end if
             end do
 
             ! otherwise, just copy element-to-element
          else
             do i = 1, s_target_list%size()
                call space_interp%map(s_target_list%x(i), this%s(i)%x, &
                     this%nelv, xh)
             end do
          end if ! present s_index_list
       end if ! present s_tgt
 
       call space_interp%free
       call prev_xh%free
 
    end if ! if interpolate
 
    nullify(dof)
    nullify(xh)
    nullify(msh)
 
  end subroutine fld_file_data_import_fields
 
  subroutine fld_file_data_init(this, nelv, offset_el)
    class(fld_file_data_t), intent(inout) :: this
    integer, intent(in), optional :: nelv, offset_el
 
    call this%free()
    if (present(nelv)) this%nelv = nelv
    if (present(offset_el)) this%offset_el = offset_el
 
  end subroutine fld_file_data_init
 
  function fld_file_data_size(this) result(i)
    class(fld_file_data_t) :: this
    integer :: i
    i = 0
    if (this%u%size() .gt. 0) i = i + 1
    if (this%v%size() .gt. 0) i = i + 1
    if (this%w%size() .gt. 0) i = i + 1
    if (this%p%size() .gt. 0) i = i + 1
    if (this%t%size() .gt. 0) i = i + 1
    i = i + this%n_scalars
 
  end function fld_file_data_size
 
  subroutine fld_file_data_init_same(this, fld_file, n)
    class(fld_file_data_t), target, intent(inout) :: this
    class(fld_file_data_t), target, intent(in) :: fld_file
    integer, intent(in) :: n
    integer :: i, j
 
    if (fld_file%u%size() .gt. 0) then
       call this%u%init(n)
    end if
    if (fld_file%v%size() .gt. 0) then
       call this%v%init(n)
    end if
    if (fld_file%w%size() .gt. 0) then
       call this%w%init(n)
    end if
    if (fld_file%p%size() .gt. 0) then
       call this%p%init(n)
    end if
    if (fld_file%t%size() .gt. 0) then
       call this%t%init(n)
    end if
    this%n_scalars = fld_file%n_scalars
    allocate(this%s(fld_file%n_scalars))
    do j = 1, fld_file%n_scalars
       call this%s(j)%init(n)
    end do
 
  end subroutine fld_file_data_init_same
 
  subroutine fld_file_data_init_n_fields(this, n_fields, n)
    class(fld_file_data_t), target, intent(inout) :: this
    integer, intent(in) :: n, n_fields
    integer :: i, j
 
 
    if (n_fields .gt. 0) then
       call this%u%init(n)
    end if
    if (n_fields .gt. 1) then
       call this%v%init(n)
    end if
    if (n_fields .gt. 2) then
       call this%w%init(n)
    end if
    if (n_fields .gt. 3) then
       call this%p%init(n)
    end if
    if (n_fields .gt. 4) then
       call this%t%init(n)
    end if
    if (n_fields .gt. 5) then
       this%n_scalars = n_fields-5
       allocate(this%s(this%n_scalars))
       do j = 1, this%n_scalars
          call this%s(j)%init(n)
       end do
    end if
 
  end subroutine fld_file_data_init_n_fields
 
  subroutine fld_file_data_get_list(this, ptr_list, n)
    class(fld_file_data_t), target, intent(in) :: this
    integer, intent(in) :: n
    integer :: i, j
    type(vector_ptr_t), intent(inout) :: ptr_list(n)
    i = 1
    if (this%u%size() .gt. 0) then
       ptr_list(i)%ptr => this%u
       i = i + 1
    end if
    if (this%v%size() .gt. 0) then
       ptr_list(i)%ptr => this%v
       i = i + 1
    end if
    if (this%w%size() .gt. 0) then
       ptr_list(i)%ptr => this%w
       i = i + 1
    end if
    if (this%p%size() .gt. 0) then
       ptr_list(i)%ptr => this%p
       i = i + 1
    end if
    if (this%t%size() .gt. 0) then
       ptr_list(i)%ptr => this%t
       i = i + 1
    end if
    do j = 1, this%n_scalars
       ptr_list(i)%ptr => this%s(j)
       i = i +1
    end do
 
  end subroutine fld_file_data_get_list
 
 
 
  subroutine fld_file_data_scale(this, c)
    class(fld_file_data_t), intent(inout) :: this
    real(kind=rp), intent(in) :: c
    integer :: i
 
    if (this%u%size() .gt. 0) call cmult(this%u%x, c, this%u%size())
    if (this%v%size() .gt. 0) call cmult(this%v%x, c, this%v%size())
    if (this%w%size() .gt. 0) call cmult(this%w%x, c, this%w%size())
    if (this%p%size() .gt. 0) call cmult(this%p%x, c, this%p%size())
    if (this%t%size() .gt. 0) call cmult(this%t%x, c, this%t%size())
 
    do i = 1, this%n_scalars
       if (this%s(i)%size() .gt. 0) call cmult(this%s(i)%x, c, this%s(i)%size())
    end do
 
  end subroutine fld_file_data_scale
 
  subroutine fld_file_data_add(this, other)
    class(fld_file_data_t), intent(inout) :: this
    class(fld_file_data_t), intent(in) :: other
    integer :: i, n
 
    if (this%u%size() .gt. 0) call add2(this%u%x, other%u%x, this%u%size())
    if (this%v%size() .gt. 0) call add2(this%v%x, other%v%x, this%v%size())
    if (this%w%size() .gt. 0) call add2(this%w%x, other%w%x, this%w%size())
    if (this%p%size() .gt. 0) call add2(this%p%x, other%p%x, this%p%size())
    if (this%t%size() .gt. 0) call add2(this%t%x, other%t%x, this%t%size())
 
    do i = 1, this%n_scalars
       if (this%s(i)%size() .gt. 0) call add2(this%s(i)%x, other%s(i)%x, &
            this%s(i)%size())
    end do
  end subroutine fld_file_data_add
 
  subroutine fld_file_data_free(this)
    class(fld_file_data_t), intent(inout) :: this
    integer :: i
    call this%x%free()
    call this%y%free()
    call this%z%free()
    call this%u%free()
    call this%v%free()
    call this%w%free()
    call this%p%free()
    call this%t%free()
    if (allocated(this%s)) then
       do i = 1, this%n_scalars
          call this%s(i)%free()
       end do
       deallocate(this%s)
    end if
    this%n_scalars = 0
    this%time = 0.0
    this%glb_nelv = 0
    this%nelv = 0
    this%offset_el = 0
    this%lx = 0
    this%ly = 0
    this%lz = 0
    this%t_counter = 0
    this%meta_nsamples = 0
    this%meta_start_counter = 0
    if (allocated(this%idx)) deallocate(this%idx)
  end subroutine fld_file_data_free
 
  subroutine fld_file_data_generate_interpolator(this, global_interp, to_dof, &
       to_msh, tolerance, padding, global_interp_subdict)
    class(fld_file_data_t), intent(in) :: this
    type(global_interpolation_t), intent(inout) :: global_interp
    type(dofmap_t), intent(in), target :: to_dof
    type(mesh_t), intent(in), target :: to_msh
    real(kind=dp), intent(in), optional :: tolerance
    real(kind=dp), intent(in), optional :: padding
    type(json_file), intent(inout), optional :: global_interp_subdict
 
    ! --- variables for interpolation
    type(space_t) :: fld_Xh
    real(kind=rp), allocatable :: x_coords(:,:,:,:), y_coords(:,:,:,:), &
         z_coords(:,:,:,:)
    real(kind=rp) :: center_x, center_y, center_z
    integer :: e, i
    real(kind=rp) :: tol_ = glob_interp_tol
    ! ---
 
    type(space_t), pointer :: to_xh
    to_xh => to_dof%Xh
 
    ! Safeguard in case we didn't read mesh information
    if (.not. allocated(this%x%x) .or. &
         .not. allocated(this%y%x) .or. &
         .not. allocated(this%z%x)) call neko_error("Unable to retrieve &
    &mesh information from fld data.")
 
    ! Create a space based on the fld data
    call fld_xh%init(gll, this%lx, this%ly, this%lz)
 
    ! These are the coordinates of our current dofmap
    ! that we use for the interpolation
    allocate(x_coords(to_xh%lx, to_xh%ly, to_xh%lz, to_msh%nelv))
    allocate(y_coords(to_xh%lx, to_xh%ly, to_xh%lz, to_msh%nelv))
    allocate(z_coords(to_xh%lx, to_xh%ly, to_xh%lz, to_msh%nelv))
 
    ! JSON subdict takes precedence over dummy argument, but if neither is
    ! provided throw an error
    if (present(global_interp_subdict)) then
       call json_get_or_default(global_interp_subdict, "tolerance", tol_, &
            tol_)
    else if (present(tolerance)) then
       tol_ = tolerance
    else
       call neko_error("No tolerance provided for interpolation! " // &
            "Please provide a tolerance or a subdict with a tolerance entry.")
    end if
 
    do e = 1, to_msh%nelv
       center_x = 0d0
       center_y = 0d0
       center_z = 0d0
       do i = 1, to_xh%lxyz
          center_x = center_x + to_dof%x(i, 1, 1, e)
          center_y = center_y + to_dof%y(i, 1, 1, e)
          center_z = center_z + to_dof%z(i, 1, 1, e)
       end do
       center_x = center_x / to_xh%lxyz
       center_y = center_y / to_xh%lxyz
       center_z = center_z / to_xh%lxyz
       do i = 1, to_xh%lxyz
          x_coords(i, 1, 1, e) = to_dof%x(i, 1, 1, e) - &
               tol_ * (to_dof%x(i, 1, 1, e) - center_x)
          y_coords(i, 1, 1, e) = to_dof%y(i, 1, 1, e) - &
               tol_ * (to_dof%y(i, 1, 1, e) - center_y)
          z_coords(i, 1, 1, e) = to_dof%z(i, 1, 1, e) - &
               tol_ * (to_dof%z(i, 1, 1, e) - center_z)
       end do
    end do
 
    ! The initialization is done based on the variables created from
    ! fld data.
    if (present(global_interp_subdict)) then
       call global_interp%init(this%x%x, this%y%x, this%z%x, &
            this%gdim, this%nelv, fld_xh, &
            global_interp_subdict)
    else
       call global_interp%init(this%x%x, this%y%x, this%z%x, &
            this%gdim, this%nelv, fld_xh, &
            tol = tolerance, pad = padding)
    end if
 
    call global_interp%find_points(x_coords, y_coords, z_coords, &
         to_dof%size())
 
    deallocate(x_coords)
    deallocate(y_coords)
    deallocate(z_coords)
    call fld_xh%free()
 
  end subroutine fld_file_data_generate_interpolator
 
end module fld_file_data