cell_m.f90

!------------------------------------------------------------------------------------------------------------------------------------
! Copyright  2018 Fabrice Roy
!
! Contact: fabrice.roy@obspm.fr
!
! This file is part of pFoF.
!
! pFoF is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! pFoF is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with pFoF. If not, see <http://www.gnu.org/licenses/>.

!> @file
!! Class for Ramses AMR cell
!! @brief
!!
!! @author Fabrice Roy

!> Class for Ramses AMR cell
!------------------------------------------------------------------------------------------------------------------------------------

module cell_m

  use constants_m, only : DP, EPSILON_D, IDKIND, MPI_DP, MPI_IDKIND
  use iso_fortran_env, only : ERROR_UNIT, OUTPUT_UNIT
  use shared_variables_m

  implicit none

  private

  public :: cell_t, Init_cell_mpi_type

  type cell_t(n_non_thermal_pressures, n_passive_scalars)
     integer, len :: n_non_thermal_pressures
     integer, len :: n_passive_scalars
     integer(kind=4) :: level
     integer(kind=4) :: process_id
     integer(kind=4) :: cube_id
     integer(kind=IDKIND) :: identity
     real(kind=DP), dimension(3) :: position
     real(kind=DP), dimension(3) :: velocity_field
     real(kind=DP) :: density
     real(kind=DP) :: thermal_pressure
     real(kind=DP), dimension(n_non_thermal_pressures) :: non_thermal_pressures
     real(kind=DP), dimension(n_passive_scalars) :: passive_scalars
   contains
     procedure :: Is_eq
     procedure :: Is_sup
     procedure :: Print => Print_cell
     procedure :: Write => Write_hdf5_cell
  end type cell_t

contains

  !------------------------------------------------------------------------------------------------------------------------------------
  subroutine Init_cell_mpi_type(cell_mpi_t, cell_array) 

    use mpi

    integer, intent(out) :: cell_mpi_t
    class(cell_t(*,*)), intent(in), dimension(*) :: cell_array

    integer :: im
    integer :: mpierr
    integer(kind=MPI_ADDRESS_KIND) :: displs(12), lower_bound, extent
    integer :: blocklengths(12)=(/1,1,1,1,1,1,3,3,1,1,0,0/),         &
         types(12)=(/MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_IDKIND, MPI_DP, MPI_DP, MPI_DP, MPI_DP, MPI_DP, MPI_DP/)
    integer :: tmp_type

    call timer%Set_ref()
    call mpi_get_address(cell_array(1)%n_non_thermal_pressures, displs(1), mpierr)
    call mpi_get_address(cell_array(1)%n_passive_scalars, displs(2), mpierr)
    call mpi_get_address(cell_array(1)%level, displs(3), mpierr)
    call mpi_get_address(cell_array(1)%process_id, displs(4), mpierr)
    call mpi_get_address(cell_array(1)%cube_id, displs(5), mpierr)
    call mpi_get_address(cell_array(1)%identity, displs(6), mpierr)
    call mpi_get_address(cell_array(1)%position(1), displs(7), mpierr)
    call mpi_get_address(cell_array(1)%velocity_field(1), displs(8), mpierr)
    call mpi_get_address(cell_array(1)%density, displs(9), mpierr)
    call mpi_get_address(cell_array(1)%thermal_pressure, displs(10), mpierr)
    call mpi_get_address(cell_array(1)%non_thermal_pressures, displs(11), mpierr)
    call mpi_get_address(cell_array(1)%passive_scalars, displs(12), mpierr)
    call mpi_get_address(cell_array(1), lower_bound, mpierr)
    call mpi_get_address(cell_array(2), extent, mpierr)

    do im = 1, 12
       displs(im) = displs(im) - lower_bound
    end do
    extent = extent - lower_bound
    lower_bound = 0

    call mpi_type_create_struct(12, blocklengths, displs, types, tmp_type, mpierr)
    call mpi_type_commit(tmp_type, mpierr)
    call mpi_type_create_resized(tmp_type, lower_bound, extent, cell_mpi_t, mpierr)
    call mpi_type_free(tmp_type, mpierr)
    call mpi_type_commit(cell_mpi_t, mpierr)

    call timer%Inc_comm()
  end subroutine Init_cell_mpi_type

#ifdef DEBUG2
  write(ERROR_UNIT,'(a,i0)') 'Is_eq begins on process', mpi_process%rank
#endif

  !------------------------------------------------------------------------------------------------------------------------------------
  pure function Is_eq(this,other_cell) result (is_equal)

    class(cell_t(*,*)), intent(in) :: this
    class(cell_t(*,*)), intent(in) :: other_cell

    logical :: is_equal

#ifdef DEBUG2
    write(ERROR_UNIT,'(a,i0)') 'Is_eq begins on process', mpi_process%rank
#endif

    is_equal = .true.
    if(this%process_id /= other_cell%process_id) then
       is_equal = .false.
       return
    else if(this%cube_id /= other_cell%cube_id) then
       is_equal = .false.
       return
    else if( this%level /= other_cell%level) then
       is_equal = .false.
       return
    else if( this%position(1) - other_cell%position(1) > EPSILON_D ) then
       is_equal = .false.
       return
    else if( this%position(2) - other_cell%position(2) > EPSILON_D ) then
       is_equal = .false.
       return
    else if( this%position(3) - other_cell%position(3) > EPSILON_D ) then
       is_equal = .false.
       return
    end if

#ifdef DEBUG2
    write(ERROR_UNIT,'(a,i0)') 'Is_eq ends on process', mpi_process%rank
#endif

  end function Is_eq


  !------------------------------------------------------------------------------------------------------------------------------------
  pure function Is_sup(this,other_cell) result (is_superior)

    class(cell_t(*,*)), intent(in) :: this
    class(cell_t(*,*)), intent(in) :: other_cell

    logical :: is_superior

#ifdef DEBUG2
    write(ERROR_UNIT,'(a,i0)') 'Is_sup begins on process', mpi_process%rank
#endif

    is_superior = .false.
    ! superior means: higher process_id, higher cube_id, higher level, then higher x, then higher y, then higher z
    if( this%process_id > other_cell%process_id) then
       is_superior = .true.
       return
    else if (this%process_id < other_cell%process_id) then
       is_superior = .false.
       return
    else if (this%process_id == other_cell%process_id) then
       if ( this%cube_id > other_cell%cube_id ) then
          is_superior = .true.
          return
       else if( this%cube_id < other_cell%cube_id ) then
          is_superior = .false. 
          return
       else if ( this%cube_id == other_cell%cube_id ) then
          if( this%level > other_cell%level ) then
             is_superior = .true.
             return
          else if ( this%level < other_cell%level ) then
             is_superior = .false.
             return
          else if ( this%level == other_cell%level ) then
             if ( this%position(1) > other_cell%position(1) ) then
                is_superior = .true.
                return
             else if ( this%position(1) < other_cell%position(1) ) then
                is_superior = .false.
                return
             else if ( this%position(1) - other_cell%position(1) < EPSILON_D ) then
                if ( this%position(2) > other_cell%position(2) ) then
                   is_superior = .true.
                   return
                else if ( this%position(2) < other_cell%position(2) ) then
                   is_superior = .false.
                   return
                else if ( this%position(2) - other_cell%position(2) < EPSILON_D) then
                   if ( this%position(3) > other_cell%position(3) ) then
                      is_superior = .true.
                      return
                   else if ( this%position(3) < other_cell%position(3) ) then
                      is_superior = .false.
                      return
                   end if
                end if
             end if
          end if
       end if
    end if

#ifdef DEBUG2
    write(ERROR_UNIT,'(a,i0)') 'Is_sup ends on process', mpi_process%rank
#endif

  end function Is_sup


  !------------------------------------------------------------------------------------------------------------------------------------
  subroutine Print_cell(this, unit_number)

    class(cell_t(*,*)), intent(in) :: this
    integer, intent(in) :: unit_number

#ifdef DEBUG
    write(ERROR_UNIT,'(a,i0)') 'Print_cell begins on process', mpi_process%rank
#endif
    call timer%Set_ref()

    write(unit_number,'(a,i0)') 'cell:', this%identity
    write(unit_number,'(a,i0)') '- level:', this%level
    write(unit_number,'(a,i0)') '- process_id:', this%process_id
    write(unit_number,'(a,i0)') '- cube_id:', this%cube_id
    write(unit_number,'(a,f10.8,a,f10.8,a,f10.8,a)') '- position: (', this%position(1),',',this%position(2),',',this%position(3),')'
    write(unit_number,'(a,f10.8,a,f10.8,a,f10.8,a)') '- velocity_field: (', this%velocity_field(1),',',this%velocity_field(2),',',this%velocity_field(3),')'
    write(unit_number,*) '- thermal_pressure:', this%thermal_pressure
    write(unit_number,'(a,i0)') '- n_non_thermal_pressures:', this%n_non_thermal_pressures
    if(this%n_non_thermal_pressures /= 0) write(unit_number,*) '- non_thermal_pressures:', this%non_thermal_pressures(:)
    write(unit_number,'(a,i0)') '- n_passive_scalars:', this%n_passive_scalars
    if(this%n_passive_scalars /= 0) write(unit_number,*) '- passive_scalars:', this%passive_scalars(:)

    call timer%Inc_out()
#ifdef DEBUG
    write(ERROR_UNIT,'(a,i0)') 'Print_cell ends on process', mpi_process%rank
#endif

  end subroutine Print_cell


  !------------------------------------------------------------------------------------------------------------------------------------
  subroutine Write_hdf5_cell(this, hdf5_id)

    use char_utils_m, only : Int_to_char12, Int_to_char5
    use fortran_hdf5_constants_m
    use fortran_hdf5_manage_groups_m
    use fortran_hdf5_write_attribute_m
    use fortran_hdf5_write_data_m
    use hdf5

    class(cell_t(*,*)),intent(in) :: this
    integer(kind=hid_t),intent(in) :: hdf5_id

    character(len=H5_STR_LEN) :: attr_name
    character(len=H5_STR_LEN) :: groupname
    integer(kind=hid_t) :: cell_gr_id
    integer :: idim

#ifdef DEBUG
    write(ERROR_UNIT,'(a,i0)') 'Write_hdf5_cell begins on process', mpi_process%rank
#endif
    call timer%Set_ref()

    attr_name = 'identity'
    call Hdf5_Write_data(hdf5_id, attr_name, this%identity)
    attr_name = 'level'
    call Hdf5_Write_attr(hdf5_id, attr_name, this%level)
    attr_name = 'position'
    call Hdf5_Write_attr(hdf5_id, attr_name, 3, this%position)
    attr_name = 'velocity field'
    call Hdf5_Write_attr(hdf5_id, attr_name, 3, this%velocity_field)
    attr_name = 'density'
    call Hdf5_Write_attr(hdf5_id, attr_name, this%density)
    attr_name = 'thermal pressure'
    call Hdf5_Write_attr(hdf5_id, attr_name, this%thermal_pressure)
    if ( this%n_non_thermal_pressures /= 0 ) then
       attr_name = 'non thermal pressures'
       call Hdf5_Write_attr(hdf5_id, attr_name, this%n_non_thermal_pressures, this%non_thermal_pressures)
    end if
    if( this%n_passive_scalars /= 0 ) then 
       attr_name = 'passive scalars'
       call Hdf5_Write_attr(hdf5_id, attr_name, this%n_passive_scalars, this%passive_scalars)
    end if

    call timer%Inc_out()
#ifdef DEBUG
    write(ERROR_UNIT,'(a,i0)') 'Write_hdf5_cell ends on process', mpi_process%rank
#endif

  end subroutine Write_hdf5_cell


end module cell_m