MPI-AMRVAC  3.0
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_convert.t
Go to the documentation of this file.
1 module mod_convert
2 
3  use mpi
4  use mod_variables, only: max_nw
5 
6  implicit none
7  public
8 
9  abstract interface
10 
11  function sub_convert_vars(ixI^L, ixO^L, w, x, nwc) result(wnew)
13  integer, intent(in) :: ixi^l,ixo^l, nwc
14  double precision, intent(in) :: w(ixi^s, 1:nw)
15  double precision, intent(in) :: x(ixi^s,1:ndim)
16  double precision :: wnew(ixo^s, 1:nwc)
17  end function sub_convert_vars
18 
19 
20  end interface
22  procedure(sub_convert_vars), pointer, nopass :: phys_convert_vars
23  integer :: nwc
24  character(len=40) :: file_suffix
25  character(len=40) :: dataset_names(max_nw)
26  type(convert_vars_method), pointer :: next
27  end type convert_vars_method
29 contains
30 
31 
32  subroutine init_convert()
33  !initialize the head of the linked list of convert methods
35  end subroutine init_convert
36 
37 
38 
39  function get_names_from_string(aux_variable_names,nwc) result(names)
41  character(len=*),intent(in):: aux_variable_names
42  integer, intent(in) :: nwc
43  character(len=name_len) :: names(1:nwc)
44 
45  integer:: space_position,iw
46  character(len=name_len):: wname
47  character(len=std_len):: scanstring
48 
49 
50  ! copied from subroutine getheadernames in calculate_xw
51  scanstring=trim(adjustl(aux_variable_names))
52  space_position=index(scanstring,' ')
53  do iw=1,nwc
54  do while (space_position==1)
55  scanstring=scanstring(2:)
56  space_position=index(scanstring,' ')
57  enddo
58  wname=scanstring(:space_position-1)
59  scanstring=scanstring(space_position+1:)
60  space_position=index(scanstring,' ')
61 
62  names(iw)=trim(adjustl(wname))
63  enddo
64  end function get_names_from_string
65 
66 
67  ! shortcut
68  subroutine add_convert_method2(phys_convert_vars, nwc, aux_variable_names, file_suffix)
70  integer, intent(in) :: nwc
71  character(len=*),intent(in):: aux_variable_names
72  character(len=*), intent(in) :: file_suffix
73 
74  interface
75  function phys_convert_vars(ixI^L, ixO^L, w, x, nwc) result(wnew)
77  integer, intent(in) :: ixI^L, ixO^L, nwc
78  double precision, intent(in) :: w(ixI^S, 1:nw)
79  double precision, intent(in) :: x(ixI^S,1:ndim)
80  double precision :: wnew(ixO^S, 1:nwc)
81  end function phys_convert_vars
82  end interface
83 
84  call add_convert_method(phys_convert_vars, nwc, get_names_from_string(aux_variable_names,nwc), file_suffix)
85 
86  end subroutine add_convert_method2
87 
88  subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
89  integer, intent(in) :: nwc
90  character(len=*), intent(in) :: dataset_names(:)
91  character(len=*), intent(in) :: file_suffix
92 
93  interface
94  function phys_convert_vars(ixI^L, ixO^L, w, x, nwc) result(wnew)
96  integer, intent(in) :: ixI^L, ixO^L, nwc
97  double precision, intent(in) :: w(ixI^S, 1:nw)
98  double precision, intent(in) :: x(ixI^S,1:ndim)
99  double precision :: wnew(ixO^S, 1:nwc)
100  end function phys_convert_vars
101  end interface
102 
103  type(convert_vars_method), pointer :: temp
104 
105  if(nwc .gt. max_nw) then
106  call mpistop("INCREASE max_nw ")
107  endif
108 
109  allocate(temp)
110  temp%phys_convert_vars => phys_convert_vars
111  temp%nwc = nwc
112  temp%file_suffix = file_suffix
113  temp%dataset_names(1:nwc) = dataset_names(1:nwc)
114  temp%next => head_convert_vars_methods
116  end subroutine add_convert_method
117 
118 
119  subroutine convert_all()
120  type(convert_vars_method), pointer :: temp
122  do while(associated(temp))
123  call convert_dat_generic(temp%nwc, temp%dataset_names, temp%file_suffix, temp%phys_convert_vars)
124  temp=>temp%next
125  end do
126  end subroutine convert_all
127 
128  ! Copied from subroutine write_snapshot in amrvacio/mod_input_output
129  ! the staggered values are not saved in this subroutine!
130  subroutine convert_dat_generic(nwc, dataset_names, file_suffix, convert_vars)
131 
132  use mod_forest
135 
136  integer :: file_handle, igrid, Morton_no, iwrite
137  integer :: ipe, ix_buffer(2*ndim+1), n_values
138  integer :: ixI^L, ixO^L, n_ghost(2*ndim)
139  integer :: ixOs^L,n_values_stagger
140  integer :: iorecvstatus(MPI_STATUS_SIZE)
141  integer :: ioastatus(MPI_STATUS_SIZE)
142  integer :: igrecvstatus(MPI_STATUS_SIZE)
143  integer :: istatus(MPI_STATUS_SIZE)
144  type(tree_node), pointer :: pnode
145  integer(kind=MPI_OFFSET_KIND) :: offset_tree_info
146  integer(kind=MPI_OFFSET_KIND) :: offset_block_data
147  integer(kind=MPI_OFFSET_KIND) :: offset_offsets
148  double precision, allocatable :: w_buffer(:)
149  double precision, allocatable:: converted_vars(:^D&,:)
150 
151  integer :: idim, itag
152  integer, allocatable :: block_ig(:, :)
153  integer, allocatable :: block_lvl(:)
154  integer(kind=MPI_OFFSET_KIND), allocatable :: block_offset(:)
155 
156  integer, intent(in) :: nwc
157  character(len=*), intent(in) :: dataset_names(:)
158  character(len=*), intent(in) :: file_suffix
159  interface
160 
161  function convert_vars(ixI^L, ixO^L,w,x,nwc) result(wres)
163  integer, intent(in) :: ixI^L, ixO^L, nwc
164  double precision, intent(in) :: x(ixI^S,1:ndim)
165  double precision, intent(in) :: w(ixI^S,1:nw)
166  double precision :: wres(ixO^S,1:nwc)
167  end function convert_vars
168 
169  end interface
170 
171 
172  call mpi_barrier(icomm, ierrmpi)
173 
174  ! Allocate send/receive buffer
175  n_values = count_ix(ixg^ll) * nw
176  allocate(w_buffer(n_values))
177 
178  ! Allocate arrays with information about grid blocks
179  allocate(block_ig(ndim, nleafs))
180  allocate(block_lvl(nleafs))
181  allocate(block_offset(nleafs+1))
182 
183  ! master processor
184  if (mype==0) then
185  call create_output_file(file_handle, snapshotnext, ".dat", file_suffix)
186 
187  ! Don't know offsets yet, we will write header again later
188  offset_tree_info = -1
189  offset_block_data = -1
190  call snapshot_write_header1(file_handle, offset_tree_info, &
191  offset_block_data, dataset_names, nwc)
192 
193  call mpi_file_get_position(file_handle, offset_tree_info, ierrmpi)
194 
195  call write_forest(file_handle)
196 
197  ! Collect information about the spatial index (ig^D) and refinement level
198  ! of leaves
199  do morton_no = morton_start(0), morton_stop(npe-1)
200  igrid = sfc(1, morton_no)
201  ipe = sfc(2, morton_no)
202  pnode => igrid_to_node(igrid, ipe)%node
203 
204  block_ig(:, morton_no) = [ pnode%ig^d ]
205  block_lvl(morton_no) = pnode%level
206  block_offset(morton_no) = 0 ! Will be determined later
207  end do
208 
209  call mpi_file_write(file_handle, block_lvl, size(block_lvl), &
210  mpi_integer, istatus, ierrmpi)
211 
212  call mpi_file_write(file_handle, block_ig, size(block_ig), &
213  mpi_integer, istatus, ierrmpi)
214 
215  ! Block offsets are currently unknown, but will be overwritten later
216  call mpi_file_get_position(file_handle, offset_offsets, ierrmpi)
217  call mpi_file_write(file_handle, block_offset(1:nleafs), nleafs, &
218  mpi_offset, istatus, ierrmpi)
219 
220  call mpi_file_get_position(file_handle, offset_block_data, ierrmpi)
221 
222  ! Check whether data was written as expected
223  if (offset_block_data - offset_tree_info /= &
224  (nleafs + nparents) * size_logical + &
225  nleafs * ((1+ndim) * size_int + 2 * size_int)) then
226  if (mype == 0) then
227  print *, "Warning: MPI_OFFSET type /= 8 bytes"
228  print *, "This *could* cause problems when reading .dat files"
229  end if
230  end if
231 
232  block_offset(1) = offset_block_data
233  iwrite = 0
234  end if
235 
236  do morton_no=morton_start(mype), morton_stop(mype)
237  igrid = sfc_to_igrid(morton_no)
238  itag = morton_no
239  block=>ps(igrid)
240  ! this might be used in convert function,
241  ! it was not used when the output is already computed vars (write_snapshot)
242  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
243 
244  ! start copied from block_shape_io,
245  ! because nwc is needed as parameter
246  ! TODO check if this will be used elsewehere and put it in separate subroutine
247  n_ghost(:) = 0
248 
249  if(save_physical_boundary) then
250  do idim=1,ndim
251  ! Include ghost cells on lower boundary
252  if(ps(igrid)%is_physical_boundary(2*idim-1)) n_ghost(idim)=nghostcells
253  ! Include ghost cells on upper boundary
254  if(ps(igrid)%is_physical_boundary(2*idim)) n_ghost(ndim+idim)=nghostcells
255  end do
256  end if
257 
258  {ixomin^d = ixmlo^d - n_ghost(^d)\}
259  {ixomax^d = ixmhi^d + n_ghost(ndim+^d)\}
260 
261  n_values = count_ix(ixo^l) * nwc
262  ! end copied from block_shape_io
263 
264  {iximin^d = ixglo^d\}
265  {iximax^d = ixghi^d\}
266 
267  w_buffer(1:n_values) = pack(convert_vars(ixi^l, ixo^l,ps(igrid)%w(ixi^s, 1:nw), ps(igrid)%x(ixi^s, 1:ndim),nwc), .true.)
268 
269  ix_buffer(1) = n_values
270  ix_buffer(2:) = n_ghost
271 
272  if (mype /= 0) then
273  call mpi_send(ix_buffer, 2*ndim+1, &
274  mpi_integer, 0, itag, icomm, ierrmpi)
275  call mpi_send(w_buffer, n_values, &
276  mpi_double_precision, 0, itag, icomm, ierrmpi)
277  else
278  iwrite = iwrite+1
279  call mpi_file_write(file_handle, ix_buffer(2:), &
280  2*ndim, mpi_integer, istatus, ierrmpi)
281  call mpi_file_write(file_handle, w_buffer, &
282  n_values, mpi_double_precision, istatus, ierrmpi)
283 
284  ! Set offset of next block
285  block_offset(iwrite+1) = block_offset(iwrite) + &
286  int(n_values, mpi_offset_kind) * size_double + &
287  2 * ndim * size_int
288  end if
289  end do
290 
291  ! Write data communicated from other processors
292  if (mype == 0) then
293  do ipe = 1, npe-1
294  do morton_no=morton_start(ipe), morton_stop(ipe)
295  iwrite=iwrite+1
296  itag=morton_no
297 
298  call mpi_recv(ix_buffer, 2*ndim+1, mpi_integer, ipe, itag, icomm,&
299  igrecvstatus, ierrmpi)
300  n_values = ix_buffer(1)
301 
302  call mpi_recv(w_buffer, n_values, mpi_double_precision,&
303  ipe, itag, icomm, iorecvstatus, ierrmpi)
304 
305  call mpi_file_write(file_handle, ix_buffer(2:), &
306  2*ndim, mpi_integer, istatus, ierrmpi)
307  call mpi_file_write(file_handle, w_buffer, &
308  n_values, mpi_double_precision, istatus, ierrmpi)
309 
310  ! Set offset of next block
311  block_offset(iwrite+1) = block_offset(iwrite) + &
312  int(n_values, mpi_offset_kind) * size_double + &
313  2 * ndim * size_int
314  end do
315  end do
316 
317  ! Write block offsets (now we know them)
318  call mpi_file_seek(file_handle, offset_offsets, mpi_seek_set, ierrmpi)
319  call mpi_file_write(file_handle, block_offset(1:nleafs), nleafs, &
320  mpi_offset, istatus, ierrmpi)
321 
322  ! Write header again, now with correct offsets
323  call mpi_file_seek(file_handle, 0_mpi_offset_kind, mpi_seek_set, ierrmpi)
324  call snapshot_write_header1(file_handle, offset_tree_info, &
325  offset_block_data, dataset_names, nwc)
326 
327  call mpi_file_close(file_handle, ierrmpi)
328  end if
329 
330  call mpi_barrier(icomm, ierrmpi)
331 
332  end subroutine convert_dat_generic
333 
334 end module mod_convert
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:194
subroutine write_forest(file_handle)
Definition: forest.t:283
subroutine convert_dat_generic(nwc, dataset_names, file_suffix, convert_vars)
Definition: mod_convert.t:131
subroutine add_convert_method2(phys_convert_vars, nwc, aux_variable_names, file_suffix)
Definition: mod_convert.t:69
type(convert_vars_method), pointer head_convert_vars_methods
Definition: mod_convert.t:28
character(len=name_len) function, dimension(1:nwc) get_names_from_string(aux_variable_names, nwc)
Definition: mod_convert.t:40
subroutine init_convert()
Definition: mod_convert.t:33
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Definition: mod_convert.t:89
subroutine convert_all()
Definition: mod_convert.t:120
Module with basic grid data structures.
Definition: mod_forest.t:2
integer, dimension(:), allocatable, save sfc_to_igrid
Go from a Morton number to an igrid index (for a single processor)
Definition: mod_forest.t:53
integer, dimension(:), allocatable, save morton_start
First Morton number per processor.
Definition: mod_forest.t:62
integer, save nleafs
Number of leaf block.
Definition: mod_forest.t:76
integer, dimension(:), allocatable, save morton_stop
Last Morton number per processor.
Definition: mod_forest.t:65
integer, dimension(:,:), allocatable, save sfc
Array to go from a Morton number to an igrid and processor index. Sfc(1:3, MN) contains [igrid,...
Definition: mod_forest.t:43
integer, save nparents
Number of parent blocks.
Definition: mod_forest.t:73
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
Definition: mod_forest.t:32
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
integer ixghi
Upper index of grid block arrays.
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical save_physical_boundary
True for save physical boundary cells in dat files.
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
integer ierrmpi
A global MPI error return code.
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
integer npe
The number of MPI tasks.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer nghostcells
Number of ghost cells surrounding a grid.
double precision, dimension(ndim) dxlevel
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Module for reading input and writing output.
pure integer function, public count_ix(ixOL)
Compute number of elements in index range.
subroutine, public create_output_file(fh, ix, extension, suffix)
Standard method for creating a new output file.
subroutine, public block_shape_io(igrid, n_ghost, ixOL, n_values)
Determine the shape of a block for output (whether to include ghost cells, and on which sides)
subroutine, public snapshot_write_header1(fh, offset_tree, offset_block, dataset_names, nw_vars)
Write header for a snapshot, generalize cons_wnames and nw.
integer, parameter max_nw
Maximum number of variables.
Definition: mod_variables.t:41
The data structure that contains information about a tree node/grid block.
Definition: mod_forest.t:11