MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
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 
40 
41  ! shortcut
42  subroutine add_convert_method2(phys_convert_vars, nwc, aux_variable_names, file_suffix)
45  integer, intent(in) :: nwc
46  character(len=*),intent(in):: aux_variable_names
47  character(len=*), intent(in) :: file_suffix
48 
49  interface
50  function phys_convert_vars(ixI^L, ixO^L, w, x, nwc) result(wnew)
52  integer, intent(in) :: ixI^L, ixO^L, nwc
53  double precision, intent(in) :: w(ixI^S, 1:nw)
54  double precision, intent(in) :: x(ixI^S,1:ndim)
55  double precision :: wnew(ixO^S, 1:nwc)
56  end function phys_convert_vars
57  end interface
58 
59  call add_convert_method(phys_convert_vars, nwc, get_names_from_string(aux_variable_names,nwc), file_suffix)
60 
61  end subroutine add_convert_method2
62 
63  subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
64  use mod_comm_lib, only: mpistop
65  integer, intent(in) :: nwc
66  character(len=*), intent(in) :: dataset_names(:)
67  character(len=*), intent(in) :: file_suffix
68 
69  interface
70  function phys_convert_vars(ixI^L, ixO^L, w, x, nwc) result(wnew)
72  integer, intent(in) :: ixI^L, ixO^L, nwc
73  double precision, intent(in) :: w(ixI^S, 1:nw)
74  double precision, intent(in) :: x(ixI^S,1:ndim)
75  double precision :: wnew(ixO^S, 1:nwc)
76  end function phys_convert_vars
77  end interface
78 
79  type(convert_vars_method), pointer :: temp
80 
81  if(nwc .gt. max_nw) then
82  call mpistop("INCREASE max_nw ")
83  endif
84 
85  allocate(temp)
86  temp%phys_convert_vars => phys_convert_vars
87  temp%nwc = nwc
88  temp%file_suffix = file_suffix
89  temp%dataset_names(1:nwc) = dataset_names(1:nwc)
90  temp%next => head_convert_vars_methods
92  end subroutine add_convert_method
93 
94 
95  subroutine convert_all()
96  type(convert_vars_method), pointer :: temp
98  do while(associated(temp))
99  call convert_dat_generic(temp%nwc, temp%dataset_names, temp%file_suffix, temp%phys_convert_vars)
100  temp=>temp%next
101  end do
102  end subroutine convert_all
103 
104  ! Copied from subroutine write_snapshot in amrvacio/mod_input_output
105  ! the staggered values are not saved in this subroutine!
106  subroutine convert_dat_generic(nwc, dataset_names, file_suffix, convert_vars)
107 
108  use mod_forest
112 
113  integer :: file_handle, igrid, Morton_no, iwrite
114  integer :: ipe, ix_buffer(2*ndim+1), n_values
115  integer :: ixI^L, ixO^L, n_ghost(2*ndim)
116  integer :: ixOs^L,n_values_stagger
117  integer :: iorecvstatus(MPI_STATUS_SIZE)
118  integer :: ioastatus(MPI_STATUS_SIZE)
119  integer :: igrecvstatus(MPI_STATUS_SIZE)
120  integer :: istatus(MPI_STATUS_SIZE)
121  type(tree_node), pointer :: pnode
122  integer(kind=MPI_OFFSET_KIND) :: offset_tree_info
123  integer(kind=MPI_OFFSET_KIND) :: offset_block_data
124  integer(kind=MPI_OFFSET_KIND) :: offset_offsets
125  double precision, allocatable :: w_buffer(:)
126  double precision, allocatable:: converted_vars(:^D&,:)
127 
128  integer :: idim, itag
129  integer, allocatable :: block_ig(:, :)
130  integer, allocatable :: block_lvl(:)
131  integer(kind=MPI_OFFSET_KIND), allocatable :: block_offset(:)
132 
133  integer, intent(in) :: nwc
134  character(len=*), intent(in) :: dataset_names(:)
135  character(len=*), intent(in) :: file_suffix
136  interface
137 
138  function convert_vars(ixI^L, ixO^L,w,x,nwc) result(wres)
140  integer, intent(in) :: ixI^L, ixO^L, nwc
141  double precision, intent(in) :: x(ixI^S,1:ndim)
142  double precision, intent(in) :: w(ixI^S,1:nw)
143  double precision :: wres(ixO^S,1:nwc)
144  end function convert_vars
145 
146  end interface
147 
148 
149  call mpi_barrier(icomm, ierrmpi)
150 
151  ! Allocate send/receive buffer
152  n_values = count_ix(ixg^ll) * nw
153  allocate(w_buffer(n_values))
154 
155  ! Allocate arrays with information about grid blocks
156  allocate(block_ig(ndim, nleafs))
157  allocate(block_lvl(nleafs))
158  allocate(block_offset(nleafs+1))
159 
160  ! master processor
161  if (mype==0) then
162  call create_output_file(file_handle, snapshotnext, ".dat", file_suffix)
163 
164  ! Don't know offsets yet, we will write header again later
165  offset_tree_info = -1
166  offset_block_data = -1
167  call snapshot_write_header1(file_handle, offset_tree_info, &
168  offset_block_data, dataset_names, nwc)
169 
170  call mpi_file_get_position(file_handle, offset_tree_info, ierrmpi)
171 
172  call write_forest(file_handle)
173 
174  ! Collect information about the spatial index (ig^D) and refinement level
175  ! of leaves
176  do morton_no = morton_start(0), morton_stop(npe-1)
177  igrid = sfc(1, morton_no)
178  ipe = sfc(2, morton_no)
179  pnode => igrid_to_node(igrid, ipe)%node
180 
181  block_ig(:, morton_no) = [ pnode%ig^d ]
182  block_lvl(morton_no) = pnode%level
183  block_offset(morton_no) = 0 ! Will be determined later
184  end do
185 
186  call mpi_file_write(file_handle, block_lvl, size(block_lvl), &
187  mpi_integer, istatus, ierrmpi)
188 
189  call mpi_file_write(file_handle, block_ig, size(block_ig), &
190  mpi_integer, istatus, ierrmpi)
191 
192  ! Block offsets are currently unknown, but will be overwritten later
193  call mpi_file_get_position(file_handle, offset_offsets, ierrmpi)
194  call mpi_file_write(file_handle, block_offset(1:nleafs), nleafs, &
195  mpi_offset, istatus, ierrmpi)
196 
197  call mpi_file_get_position(file_handle, offset_block_data, ierrmpi)
198 
199  ! Check whether data was written as expected
200  if (offset_block_data - offset_tree_info /= &
201  (nleafs + nparents) * size_logical + &
202  nleafs * ((1+ndim) * size_int + 2 * size_int)) then
203  if (mype == 0) then
204  print *, "Warning: MPI_OFFSET type /= 8 bytes"
205  print *, "This *could* cause problems when reading .dat files"
206  end if
207  end if
208 
209  block_offset(1) = offset_block_data
210  iwrite = 0
211  end if
212 
213  do morton_no=morton_start(mype), morton_stop(mype)
214  igrid = sfc_to_igrid(morton_no)
215  itag = morton_no
216  block=>ps(igrid)
217  ! this might be used in convert function,
218  ! it was not used when the output is already computed vars (write_snapshot)
219  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
220 
221  ! start copied from block_shape_io,
222  ! because nwc is needed as parameter
223  ! TODO check if this will be used elsewehere and put it in separate subroutine
224  n_ghost(:) = 0
225 
226  if(save_physical_boundary) then
227  do idim=1,ndim
228  ! Include ghost cells on lower boundary
229  if(ps(igrid)%is_physical_boundary(2*idim-1)) n_ghost(idim)=nghostcells
230  ! Include ghost cells on upper boundary
231  if(ps(igrid)%is_physical_boundary(2*idim)) n_ghost(ndim+idim)=nghostcells
232  end do
233  end if
234 
235  {ixomin^d = ixmlo^d - n_ghost(^d)\}
236  {ixomax^d = ixmhi^d + n_ghost(ndim+^d)\}
237 
238  n_values = count_ix(ixo^l) * nwc
239  ! end copied from block_shape_io
240 
241  {iximin^d = ixglo^d\}
242  {iximax^d = ixghi^d\}
243 
244  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.)
245 
246  ix_buffer(1) = n_values
247  ix_buffer(2:) = n_ghost
248 
249  if (mype /= 0) then
250  call mpi_send(ix_buffer, 2*ndim+1, &
251  mpi_integer, 0, itag, icomm, ierrmpi)
252  call mpi_send(w_buffer, n_values, &
253  mpi_double_precision, 0, itag, icomm, ierrmpi)
254  else
255  iwrite = iwrite+1
256  call mpi_file_write(file_handle, ix_buffer(2:), &
257  2*ndim, mpi_integer, istatus, ierrmpi)
258  call mpi_file_write(file_handle, w_buffer, &
259  n_values, mpi_double_precision, istatus, ierrmpi)
260 
261  ! Set offset of next block
262  block_offset(iwrite+1) = block_offset(iwrite) + &
263  int(n_values, mpi_offset_kind) * size_double + &
264  2 * ndim * size_int
265  end if
266  end do
267 
268  ! Write data communicated from other processors
269  if (mype == 0) then
270  do ipe = 1, npe-1
271  do morton_no=morton_start(ipe), morton_stop(ipe)
272  iwrite=iwrite+1
273  itag=morton_no
274 
275  call mpi_recv(ix_buffer, 2*ndim+1, mpi_integer, ipe, itag, icomm,&
276  igrecvstatus, ierrmpi)
277  n_values = ix_buffer(1)
278 
279  call mpi_recv(w_buffer, n_values, mpi_double_precision,&
280  ipe, itag, icomm, iorecvstatus, ierrmpi)
281 
282  call mpi_file_write(file_handle, ix_buffer(2:), &
283  2*ndim, mpi_integer, istatus, ierrmpi)
284  call mpi_file_write(file_handle, w_buffer, &
285  n_values, mpi_double_precision, istatus, ierrmpi)
286 
287  ! Set offset of next block
288  block_offset(iwrite+1) = block_offset(iwrite) + &
289  int(n_values, mpi_offset_kind) * size_double + &
290  2 * ndim * size_int
291  end do
292  end do
293 
294  ! Write block offsets (now we know them)
295  call mpi_file_seek(file_handle, offset_offsets, mpi_seek_set, ierrmpi)
296  call mpi_file_write(file_handle, block_offset(1:nleafs), nleafs, &
297  mpi_offset, istatus, ierrmpi)
298 
299  ! Write header again, now with correct offsets
300  call mpi_file_seek(file_handle, 0_mpi_offset_kind, mpi_seek_set, ierrmpi)
301  call snapshot_write_header1(file_handle, offset_tree_info, &
302  offset_block_data, dataset_names, nwc)
303 
304  call mpi_file_close(file_handle, ierrmpi)
305  end if
306 
307  call mpi_barrier(icomm, ierrmpi)
308 
309  end subroutine convert_dat_generic
310 
311 end module mod_convert
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
subroutine convert_dat_generic(nwc, dataset_names, file_suffix, convert_vars)
Definition: mod_convert.t:107
subroutine add_convert_method2(phys_convert_vars, nwc, aux_variable_names, file_suffix)
Definition: mod_convert.t:43
type(convert_vars_method), pointer head_convert_vars_methods
Definition: mod_convert.t:28
subroutine init_convert()
Definition: mod_convert.t:33
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Definition: mod_convert.t:64
subroutine convert_all()
Definition: mod_convert.t:96
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
subroutine, public write_forest(file_handle)
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)
subroutine, public create_output_file(fh, ix, extension, suffix)
Standard method for creating a new output file.
pure integer function, public count_ix(ixOL)
Compute number of elements in index range.
subroutine, public snapshot_write_header1(fh, offset_tree, offset_block, dataset_names, nw_vars)
Write header for a snapshot, generalize cons_wnames and nw.
character(len=name_len) function, dimension(1:nwc), public get_names_from_string(aux_variable_names, nwc)
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)
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