ramses_hydro_data_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 hydro data
!! @brief
!! 
!! @author Fabrice Roy

!> Class for Ramses hydro data
!------------------------------------------------------------------------------------------------------------------------------------

module ramses_hydro_data_m

  use constants_m, only : DP
  use iso_fortran_env, only : ERROR_UNIT, OUTPUT_UNIT
  use shared_variables_m

  implicit none

  private
  public :: ramses_hydro_data_t

  type ramses_hydro_data_t
     integer :: n_non_thermal_pressures
     integer :: n_passive_scalars
     real(kind=DP), allocatable, dimension(:) :: density
     real(kind=DP), allocatable, dimension(:,:) :: velocity_field
     real(kind=DP), allocatable, dimension(:,:) :: non_thermal_pressures
     real(kind=DP), allocatable, dimension(:) :: thermal_pressure
     real(kind=DP), allocatable, dimension(:,:) :: passive_scalars
   contains
     procedure :: Read => Read_ramses_hydro_data
     procedure :: Deallocate => Deallocate_ramses_hydro_data
  end type ramses_hydro_data_t

contains

  !------------------------------------------------------------------------------------------------------------------------------------
  subroutine Deallocate_ramses_hydro_data(this)

    use error_handling_m, only : Deallocate_error, ERR_MSG_LEN

    class(ramses_hydro_data_t), intent(out) :: this

    character(len=ERR_MSG_LEN) :: error_message
    integer :: alloc_stat


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

    if(allocated(this%density)) then
       write(ERROR_UNIT,*) 'Deallocate grid_index'
       deallocate(this%density, stat=alloc_stat, errmsg=error_message)
       if(alloc_stat /= 0) then
          call Deallocate_error('this%density','Deallocate_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
       end if
    end if
    if(allocated(this%velocity_field)) then
       deallocate(this%velocity_field, stat=alloc_stat, errmsg=error_message)
       if(alloc_stat /= 0) then
          call Deallocate_error('this%velocity_field','Deallocate_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
       end if
    end if
    if(allocated(this%non_thermal_pressures)) then
       deallocate(this%non_thermal_pressures, stat=alloc_stat, errmsg=error_message)
       if(alloc_stat /= 0) then
          call Deallocate_error('this%non_thermal_pressures','Deallocate_ramses_hydro_data', error_message, alloc_stat, &
               mpi_process%rank)
       end if
    end if
    if(allocated(this%thermal_pressure)) then
       deallocate(this%thermal_pressure, stat=alloc_stat, errmsg=error_message)
       if(alloc_stat /= 0) then
          call Deallocate_error('this%thermal_pressure','Deallocate_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
       end if
    end if
    if(allocated(this%passive_scalars)) then
       deallocate(this%passive_scalars, stat=alloc_stat, errmsg=error_message)
       if(alloc_stat /= 0) then
          call Deallocate_error('this%passive_scalars','Deallocate_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
       end if
    end if

    call timer%Inc_comp()
#ifdef DEBUG2
    write(ERROR_UNIT,'(a,i0)') 'Deallocate_ramses_hydro_data ends on process', mpi_process%rank
#endif

  end subroutine Deallocate_ramses_hydro_data

  !------------------------------------------------------------------------------------------------------------------------------------

  subroutine Read_ramses_hydro_data(this, ramses_hydro_header, filename, ncells)

    use error_handling_m, only : Allocate_error, IO_error, ERR_MSG_LEN
    use ramses_hydro_header_m, only : ramses_hydro_header_t

    class(ramses_hydro_data_t), intent(inout) :: this
    type(ramses_hydro_header_t), intent(inout) :: ramses_hydro_header
    character(len=*), intent(in) :: filename
    integer, intent(in) :: ncells

    integer :: ioerr
    integer :: unit_hydro
    character(len=ERR_MSG_LEN) :: error_message
    integer :: alloc_stat

    integer :: level
    integer :: twotondim
    integer :: iboundary, idim, ilevel, size, ivar
    integer :: index_begin, index_end


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

    open(newunit=unit_hydro, file=filename, status='old', action='read', form='unformatted', iostat = ioerr, iomsg = error_message)
    if( ioerr > 0 ) then
       call IO_error('open file '//filename, 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)
    end if

    read(unit_hydro, iostat=ioerr, iomsg=error_message) ramses_hydro_header%ncpu
    if(ioerr/=0) call IO_error('read ramses_hydro_header%ncpu', 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)

    read(unit_hydro, iostat=ioerr, iomsg=error_message) ramses_hydro_header%nvar
    if(ioerr/=0) call IO_error('read ramses_hydro_header%nvar', 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)

    read(unit_hydro, iostat=ioerr, iomsg=error_message) ramses_hydro_header%ndim
    if(ioerr/=0) call IO_error('read ramses_hydro_header%ndim', 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)

    read(unit_hydro, iostat=ioerr, iomsg=error_message) ramses_hydro_header%nlevel_max
    if(ioerr/=0) call IO_error('read ramses_hydro_header%nlevel_max', 'Read_ramses_hydro_file', &
         error_message, ioerr, mpi_process%rank)

    read(unit_hydro, iostat=ioerr, iomsg=error_message) ramses_hydro_header%nboundary
    if(ioerr/=0) call IO_error('read ramses_hydro_header%nboundary', 'Read_ramses_hydro_file', &
         error_message, ioerr, mpi_process%rank)

    read(unit_hydro, iostat=ioerr, iomsg=error_message) ramses_hydro_header%gamma
    if(ioerr/=0) call IO_error('read ramses_hydro_header%gamma', 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)

    ! constants def
    twotondim = 2**ramses_hydro_header%ndim

    ! allocations
    allocate(this%density(ncells), stat=alloc_stat, errmsg=error_message)
    if(alloc_stat /= 0) then
       call Allocate_error('this%density','Read_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
    end if

    allocate(this%velocity_field(ncells, ramses_hydro_header%ndim), stat=alloc_stat, errmsg=error_message)
    if(alloc_stat /= 0) then
       call Allocate_error('this%density','Read_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
    end if

    if(this%n_non_thermal_pressures/=0) then
       allocate(this%non_thermal_pressures(ncells,this%n_non_thermal_pressures), stat=alloc_stat, errmsg=error_message)
       if(alloc_stat /= 0) then
          call Allocate_error('this%density','Read_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
       end if
    end if

    allocate(this%thermal_pressure(ncells), stat=alloc_stat, errmsg=error_message)
    if(alloc_stat /= 0) then
       call Allocate_error('this%density','Read_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
    end if

    if(this%n_passive_scalars/=0) then
       allocate(this%passive_scalars(ncells,this%n_passive_scalars), stat=alloc_stat, errmsg=error_message)
       if(alloc_stat /= 0) then
          call Allocate_error('this%density','Read_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
       end if
    end if

    index_end = 0
    dolevel: do ilevel = 1, ramses_hydro_header%nlevel_max
       doboundary: do iboundary = 1, ramses_hydro_header%nboundary + ramses_hydro_header%ncpu
          read(unit_hydro, iostat=ioerr, iomsg=error_message) level
          if(ioerr/=0) call IO_error('read level', 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)

          read(unit_hydro, iostat=ioerr, iomsg=error_message) size
          if(ioerr/=0) call IO_error('read ncache', 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)

          notempty: if(size/=0) then
             index_begin = index_end + 1
             index_end = index_end + size
             dodim : do idim = 1, twotondim 
                ! read density
                read(unit_hydro, iostat=ioerr, iomsg=error_message) this%density(index_begin:index_end)
                if(ioerr/=0) call IO_error('read this%density', 'Read_ramses_hydro_file', &
                     error_message, ioerr, mpi_process%rank)
                ! read velocity
                do ivar=1, ramses_hydro_header%ndim
                   read(unit_hydro, iostat=ioerr, iomsg=error_message) this%velocity_field(index_begin:index_end,ivar)
                   if(ioerr/=0) call IO_error('read this%density', 'Read_ramses_hydro_file', &
                        error_message, ioerr, mpi_process%rank)
                end do
                ! read non thermal pressures
                if(this%n_non_thermal_pressures/=0) then
                   do ivar=1, this%n_non_thermal_pressures
                      read(unit_hydro, iostat=ioerr, iomsg=error_message) this%non_thermal_pressures(index_begin:index_end,ivar)
                      if(ioerr/=0) call IO_error('read this%non_thermal_pressures', 'Read_ramses_hydro_file', &
                           error_message, ioerr, mpi_process%rank)
                   end do
                end if
                ! read thermal pressure
                read(unit_hydro, iostat=ioerr, iomsg=error_message) this%thermal_pressure(index_begin:index_end)
                if(ioerr/=0) call IO_error('read this%thermal_pressure', 'Read_ramses_hydro_file', &
                     error_message, ioerr, mpi_process%rank)
                ! read passive scalars
                if(this%n_passive_scalars/=0) then
                   do ivar=1,this%n_passive_scalars
                      read(unit_hydro, iostat=ioerr, iomsg=error_message) this%passive_scalars(index_begin:index_end,ivar)
                      if(ioerr/=0) call IO_error('read this%passive_scalars', 'Read_ramses_hydro_file', &
                           error_message, ioerr, mpi_process%rank)
                   end do
                end if
             end do dodim
          end if notempty
       end do doboundary
    end do dolevel

    close(unit_hydro)

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

  end subroutine Read_ramses_hydro_data

end module ramses_hydro_data_m