cell_m.f90 11.3 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
20
21
22
23
24
25
26
27
28
29
30
31
!------------------------------------------------------------------------------------------------------------------------------------
! 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

32
  use constants_m, only : DP, EPSILON_D, IDKIND, MPI_DP, MPI_IDKIND
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
39
  implicit none

  private

40
  public :: cell_t, Init_cell_mpi_type
Roy Fabrice's avatar
Roy Fabrice committed
41
42
43
44
45

  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
46
     integer(kind=4) :: process_id
Roy Fabrice's avatar
Roy Fabrice committed
47
     integer(kind=4) :: cube_id
48
     integer(kind=IDKIND) :: identity
Roy Fabrice's avatar
Roy Fabrice committed
49
50
     real(kind=DP), dimension(3) :: position
     real(kind=DP), dimension(3) :: velocity_field
Roy Fabrice's avatar
Roy Fabrice committed
51
     real(kind=DP) :: density
Roy Fabrice's avatar
Roy Fabrice committed
52
53
54
55
     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
Roy Fabrice's avatar
Roy Fabrice committed
56
     procedure :: Is_eq
57
58
     procedure :: Is_sup
     procedure :: Print => Print_cell
Roy Fabrice's avatar
Roy Fabrice committed
59
     procedure :: Write => Write_hdf5_cell
Roy Fabrice's avatar
Roy Fabrice committed
60
61
62
63
  end type cell_t

contains

64
  !------------------------------------------------------------------------------------------------------------------------------------
65
  subroutine Init_cell_mpi_type(cell_mpi_t, cell_array) 
66
67
68

    use mpi

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

72
    integer :: im
73
    integer :: mpierr
Roy Fabrice's avatar
Roy Fabrice committed
74
75
76
    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/)
77
    integer :: tmp_type
78
79

    call timer%Set_ref()
80
81
82
83
    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)
Roy Fabrice's avatar
Roy Fabrice committed
84
85
86
87
88
89
90
91
    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)
92
93
94
    call mpi_get_address(cell_array(1), lower_bound, mpierr)
    call mpi_get_address(cell_array(2), extent, mpierr)

Roy Fabrice's avatar
Roy Fabrice committed
95
    do im = 1, 12
96
       displs(im) = displs(im) - lower_bound
97
98
99
100
    end do
    extent = extent - lower_bound
    lower_bound = 0

Roy Fabrice's avatar
Roy Fabrice committed
101
    call mpi_type_create_struct(12, blocklengths, displs, types, tmp_type, mpierr)
102
103
104
    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)
105
    call mpi_type_commit(cell_mpi_t, mpierr)
106
107
108
109
110
111
112

    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
Roy Fabrice's avatar
Roy Fabrice committed
113
114

  !------------------------------------------------------------------------------------------------------------------------------------
115
  pure function Is_eq(this,other_cell) result (is_equal)
Roy Fabrice's avatar
Roy Fabrice committed
116
117
118

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

Roy Fabrice's avatar
Roy Fabrice committed
120
121
    logical :: is_equal

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

Roy Fabrice's avatar
Roy Fabrice committed
126
    is_equal = .true.
127
    if(this%process_id /= other_cell%process_id) then
128
129
       is_equal = .false.
       return
Roy Fabrice's avatar
Roy Fabrice committed
130
    else if(this%cube_id /= other_cell%cube_id) then
131
132
       is_equal = .false.
       return
133
    else if( this%level /= other_cell%level) then
Roy Fabrice's avatar
Roy Fabrice committed
134
135
136
137
138
139
140
141
142
143
144
145
       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
146

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

Roy Fabrice's avatar
Roy Fabrice committed
151
152
153
  end function Is_eq


154
  !------------------------------------------------------------------------------------------------------------------------------------
155
  pure function Is_sup(this,other_cell) result (is_superior)
156
157
158

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

160
161
    logical :: is_superior

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

166
    is_superior = .false.
Roy Fabrice's avatar
Roy Fabrice committed
167
    ! superior means: higher process_id, higher cube_id, higher level, then higher x, then higher y, then higher z
168
    if( this%process_id > other_cell%process_id) then
169
170
       is_superior = .true.
       return
171
    else if (this%process_id < other_cell%process_id) then
172
173
       is_superior = .false.
       return
174
    else if (this%process_id == other_cell%process_id) then
175
       if ( this%cube_id > other_cell%cube_id ) then
176
177
          is_superior = .true.
          return
178
179
       else if( this%cube_id < other_cell%cube_id ) then
          is_superior = .false. 
180
          return
181
182
       else if ( this%cube_id == other_cell%cube_id ) then
          if( this%level > other_cell%level ) then
183
184
             is_superior = .true.
             return
185
          else if ( this%level < other_cell%level ) then
186
187
             is_superior = .false.
             return
188
189
          else if ( this%level == other_cell%level ) then
             if ( this%position(1) > other_cell%position(1) ) then
190
191
                is_superior = .true.
                return
192
             else if ( this%position(1) < other_cell%position(1) ) then
193
194
                is_superior = .false.
                return
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
             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
211
             end if
212
          end if
213
214
       end if
    end if
215
216
217
218
219

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

220
  end function Is_sup
221

222

223
224
  !------------------------------------------------------------------------------------------------------------------------------------
  subroutine Print_cell(this, unit_number)
225

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

229
230
231
#ifdef DEBUG
    write(ERROR_UNIT,'(a,i0)') 'Print_cell begins on process', mpi_process%rank
#endif
232
    call timer%Set_ref()
233
234
235

    write(unit_number,'(a,i0)') 'cell:', this%identity
    write(unit_number,'(a,i0)') '- level:', this%level
236
    write(unit_number,'(a,i0)') '- process_id:', this%process_id
Roy Fabrice's avatar
Roy Fabrice committed
237
    write(unit_number,'(a,i0)') '- cube_id:', this%cube_id
238
239
240
241
242
243
244
    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(:)
245
246

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

251
252
  end subroutine Print_cell

Roy Fabrice's avatar
Roy Fabrice committed
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274

  !------------------------------------------------------------------------------------------------------------------------------------
  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
275
    call timer%Set_ref()
Roy Fabrice's avatar
Roy Fabrice committed
276
277
278
279
280
281
282
283
284
285
286
287
288
289

    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
290
291
       attr_name = 'non thermal pressures'
       call Hdf5_Write_attr(hdf5_id, attr_name, this%n_non_thermal_pressures, this%non_thermal_pressures)
Roy Fabrice's avatar
Roy Fabrice committed
292
293
    end if
    if( this%n_passive_scalars /= 0 ) then 
294
295
       attr_name = 'passive scalars'
       call Hdf5_Write_attr(hdf5_id, attr_name, this%n_passive_scalars, this%passive_scalars)
Roy Fabrice's avatar
Roy Fabrice committed
296
    end if
297
298

    call timer%Inc_out()
Roy Fabrice's avatar
Roy Fabrice committed
299
300
301
302
303
304
305
#ifdef DEBUG
    write(ERROR_UNIT,'(a,i0)') 'Write_hdf5_cell ends on process', mpi_process%rank
#endif

  end subroutine Write_hdf5_cell


Roy Fabrice's avatar
Roy Fabrice committed
306
end module cell_m