ramses_hydro_data_m.f90 11 KB
Newer Older
Roy Fabrice's avatar
Roy Fabrice committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
!------------------------------------------------------------------------------------------------------------------------------------
! 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/>.
Roy Fabrice's avatar
Roy Fabrice committed
20

Roy Fabrice's avatar
Roy Fabrice committed
21
!> @file
Roy Fabrice's avatar
Roy Fabrice committed
22
!! Class for Ramses hydro data
Roy Fabrice's avatar
Roy Fabrice committed
23
24
25
26
!! @brief
!! 
!! @author Fabrice Roy

Roy Fabrice's avatar
Roy Fabrice committed
27
28
!> Class for Ramses hydro data
!------------------------------------------------------------------------------------------------------------------------------------
Roy Fabrice's avatar
Roy Fabrice committed
29

Roy Fabrice's avatar
Roy Fabrice committed
30
module ramses_hydro_data_m
Roy Fabrice's avatar
Roy Fabrice committed
31

Roy Fabrice's avatar
Roy Fabrice committed
32
  use constants_m, only : DP
Roy Fabrice's avatar
Roy Fabrice committed
33
  use iso_fortran_env, only : ERROR_UNIT, OUTPUT_UNIT
Roy Fabrice's avatar
Roy Fabrice committed
34
  use shared_variables_m
35

Roy Fabrice's avatar
Roy Fabrice committed
36
37
38
  implicit none

  private
Roy Fabrice's avatar
Roy Fabrice committed
39
  public :: ramses_hydro_data_t
Roy Fabrice's avatar
Roy Fabrice committed
40

Roy Fabrice's avatar
Roy Fabrice committed
41
  type ramses_hydro_data_t
42
43
44
45
46
47
48
     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
Roy Fabrice's avatar
Roy Fabrice committed
49
50
   contains
     procedure :: Read => Read_ramses_hydro_data
Roy Fabrice's avatar
Roy Fabrice committed
51
     procedure :: Deallocate => Deallocate_ramses_hydro_data
Roy Fabrice's avatar
Roy Fabrice committed
52
  end type ramses_hydro_data_t
Roy Fabrice's avatar
Roy Fabrice committed
53
54
55

contains

56
  !------------------------------------------------------------------------------------------------------------------------------------
Roy Fabrice's avatar
Roy Fabrice committed
57
  subroutine Deallocate_ramses_hydro_data(this)
Roy Fabrice's avatar
Roy Fabrice committed
58
59
60
61
62
63
64
65
66
67

    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
68
    write(ERROR_UNIT,'(a,i0)') 'Deallocate_ramses_hydro_data begins on process', mpi_process%rank
Roy Fabrice's avatar
Roy Fabrice committed
69
#endif
70
    call timer%Set_ref()
Roy Fabrice's avatar
Roy Fabrice committed
71
72
73
74
75

    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
Roy Fabrice's avatar
Roy Fabrice committed
76
          call Deallocate_error('this%density','Deallocate_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
Roy Fabrice's avatar
Roy Fabrice committed
77
78
79
80
81
       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
Roy Fabrice's avatar
Roy Fabrice committed
82
          call Deallocate_error('this%velocity_field','Deallocate_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
Roy Fabrice's avatar
Roy Fabrice committed
83
84
85
86
87
88
       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, &
Roy Fabrice's avatar
Roy Fabrice committed
89
               mpi_process%rank)
Roy Fabrice's avatar
Roy Fabrice committed
90
91
92
93
94
       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
Roy Fabrice's avatar
Roy Fabrice committed
95
          call Deallocate_error('this%thermal_pressure','Deallocate_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
Roy Fabrice's avatar
Roy Fabrice committed
96
97
98
99
100
       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
Roy Fabrice's avatar
Roy Fabrice committed
101
          call Deallocate_error('this%passive_scalars','Deallocate_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
Roy Fabrice's avatar
Roy Fabrice committed
102
103
       end if
    end if
104
105

    call timer%Inc_comp()
Roy Fabrice's avatar
Roy Fabrice committed
106
#ifdef DEBUG2
107
    write(ERROR_UNIT,'(a,i0)') 'Deallocate_ramses_hydro_data ends on process', mpi_process%rank
Roy Fabrice's avatar
Roy Fabrice committed
108
#endif
109

Roy Fabrice's avatar
Roy Fabrice committed
110
111
  end subroutine Deallocate_ramses_hydro_data

112
113
  !------------------------------------------------------------------------------------------------------------------------------------

Roy Fabrice's avatar
Roy Fabrice committed
114
  subroutine Read_ramses_hydro_data(this, ramses_hydro_header, filename, ncells)
115

116
    use error_handling_m, only : Allocate_error, IO_error, ERR_MSG_LEN
Roy Fabrice's avatar
Roy Fabrice committed
117
    use ramses_hydro_header_m, only : ramses_hydro_header_t
118

119
    class(ramses_hydro_data_t), intent(inout) :: this
Roy Fabrice's avatar
Roy Fabrice committed
120
    type(ramses_hydro_header_t), intent(inout) :: ramses_hydro_header
121
    character(len=*), intent(in) :: filename
122
    integer, intent(in) :: ncells
Roy Fabrice's avatar
Roy Fabrice committed
123
124
125

    integer :: ioerr
    integer :: unit_hydro
126
127
    character(len=ERR_MSG_LEN) :: error_message
    integer :: alloc_stat
Roy Fabrice's avatar
Roy Fabrice committed
128

Roy Fabrice's avatar
Roy Fabrice committed
129
130
    integer :: level
    integer :: twotondim
131
132
    integer :: iboundary, idim, ilevel, size, ivar
    integer :: index_begin, index_end
Roy Fabrice's avatar
Roy Fabrice committed
133
134
135


#ifdef DEBUG
136
    write(ERROR_UNIT,'(a,i0)') 'Read_ramses_hydro_data begins on process', mpi_process%rank
Roy Fabrice's avatar
Roy Fabrice committed
137
#endif
138
139
    call timer%Set_ref()

Roy Fabrice's avatar
Roy Fabrice committed
140
141
    open(newunit=unit_hydro, file=filename, status='old', action='read', form='unformatted', iostat = ioerr, iomsg = error_message)
    if( ioerr > 0 ) then
Roy Fabrice's avatar
Roy Fabrice committed
142
       call IO_error('open file '//filename, 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)
Roy Fabrice's avatar
Roy Fabrice committed
143
144
    end if

Roy Fabrice's avatar
Roy Fabrice committed
145
    read(unit_hydro, iostat=ioerr, iomsg=error_message) ramses_hydro_header%ncpu
Roy Fabrice's avatar
Roy Fabrice committed
146
    if(ioerr/=0) call IO_error('read ramses_hydro_header%ncpu', 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)
147

Roy Fabrice's avatar
Roy Fabrice committed
148
    read(unit_hydro, iostat=ioerr, iomsg=error_message) ramses_hydro_header%nvar
Roy Fabrice's avatar
Roy Fabrice committed
149
    if(ioerr/=0) call IO_error('read ramses_hydro_header%nvar', 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)
150

Roy Fabrice's avatar
Roy Fabrice committed
151
    read(unit_hydro, iostat=ioerr, iomsg=error_message) ramses_hydro_header%ndim
Roy Fabrice's avatar
Roy Fabrice committed
152
    if(ioerr/=0) call IO_error('read ramses_hydro_header%ndim', 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)
153

Roy Fabrice's avatar
Roy Fabrice committed
154
155
    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', &
Roy Fabrice's avatar
Roy Fabrice committed
156
         error_message, ioerr, mpi_process%rank)
Roy Fabrice's avatar
Roy Fabrice committed
157

Roy Fabrice's avatar
Roy Fabrice committed
158
159
    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', &
Roy Fabrice's avatar
Roy Fabrice committed
160
         error_message, ioerr, mpi_process%rank)
161

Roy Fabrice's avatar
Roy Fabrice committed
162
    read(unit_hydro, iostat=ioerr, iomsg=error_message) ramses_hydro_header%gamma
Roy Fabrice's avatar
Roy Fabrice committed
163
    if(ioerr/=0) call IO_error('read ramses_hydro_header%gamma', 'Read_ramses_hydro_file', error_message, ioerr, mpi_process%rank)
Roy Fabrice's avatar
Roy Fabrice committed
164

Roy Fabrice's avatar
Roy Fabrice committed
165
    ! constants def
Roy Fabrice's avatar
Roy Fabrice committed
166
    twotondim = 2**ramses_hydro_header%ndim
167
168
169
170

    ! allocations
    allocate(this%density(ncells), stat=alloc_stat, errmsg=error_message)
    if(alloc_stat /= 0) then
Roy Fabrice's avatar
Roy Fabrice committed
171
       call Allocate_error('this%density','Read_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
172
173
174
175
    end if

    allocate(this%velocity_field(ncells, ramses_hydro_header%ndim), stat=alloc_stat, errmsg=error_message)
    if(alloc_stat /= 0) then
Roy Fabrice's avatar
Roy Fabrice committed
176
       call Allocate_error('this%density','Read_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
177
178
179
180
181
    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
Roy Fabrice's avatar
Roy Fabrice committed
182
          call Allocate_error('this%density','Read_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
183
184
185
186
187
       end if
    end if

    allocate(this%thermal_pressure(ncells), stat=alloc_stat, errmsg=error_message)
    if(alloc_stat /= 0) then
Roy Fabrice's avatar
Roy Fabrice committed
188
       call Allocate_error('this%density','Read_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
189
190
191
192
193
    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
Roy Fabrice's avatar
Roy Fabrice committed
194
          call Allocate_error('this%density','Read_ramses_hydro_data', error_message, alloc_stat, mpi_process%rank)
195
196
       end if
    end if
197

198
    index_end = 0
Roy Fabrice's avatar
Roy Fabrice committed
199
200
    dolevel: do ilevel = 1, ramses_hydro_header%nlevel_max
       doboundary: do iboundary = 1, ramses_hydro_header%nboundary + ramses_hydro_header%ncpu
201
202
203
204
205
          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)
Roy Fabrice's avatar
Roy Fabrice committed
206

207
208
209
210
          notempty: if(size/=0) then
             index_begin = index_end + 1
             index_end = index_end + size
             dodim : do idim = 1, twotondim 
Roy Fabrice's avatar
Roy Fabrice committed
211
                ! read density
212
213
                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', &
Roy Fabrice's avatar
Roy Fabrice committed
214
                     error_message, ioerr, mpi_process%rank)
215
216
217
218
                ! 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', &
Roy Fabrice's avatar
Roy Fabrice committed
219
                        error_message, ioerr, mpi_process%rank)
220
221
222
223
224
225
                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', &
Roy Fabrice's avatar
Roy Fabrice committed
226
                           error_message, ioerr, mpi_process%rank)
227
228
229
230
231
                   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', &
Roy Fabrice's avatar
Roy Fabrice committed
232
                     error_message, ioerr, mpi_process%rank)
233
234
235
236
237
                ! 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', &
Roy Fabrice's avatar
Roy Fabrice committed
238
                           error_message, ioerr, mpi_process%rank)
239
240
                   end do
                end if
Roy Fabrice's avatar
Roy Fabrice committed
241
             end do dodim
242
          end if notempty
Roy Fabrice's avatar
Roy Fabrice committed
243
244
245
246
       end do doboundary
    end do dolevel

    close(unit_hydro)
Roy Fabrice's avatar
Roy Fabrice committed
247

248
    call timer%Inc_inp()
Roy Fabrice's avatar
Roy Fabrice committed
249
#ifdef DEBUG
250
    write(ERROR_UNIT,'(a,i0)') 'Read_ramses_hydro_data ends on process', mpi_process%rank
Roy Fabrice's avatar
Roy Fabrice committed
251
#endif
252

Roy Fabrice's avatar
Roy Fabrice committed
253
  end subroutine Read_ramses_hydro_data
254

Roy Fabrice's avatar
Roy Fabrice committed
255
end module ramses_hydro_data_m