DISPATCH
read_snapshot_mod.f90
1 !===============================================================================
2 !> IC module for reading in legacy Stagger snapshots into existing patches.
3 !>
4 !> The purpose with this separate module is to avoid having to load the entire
5 !> snapshot on each MPI rank -- there isn't enough memory to do that on the
6 !> largest experiments from the archive. In this module one can loop over either
7 !> variable index iv, or even iv and depth index, reading only one xy-slice at a
8 !> time.
9 !>
10 !> To save startup timme, separate OMP tasks are spawned for each patch, to work
11 !> in parallel on the interpolation of the ICs to the new patch mesh.
12 !>
13 !> The separate reader also saves memory because it operates before the RT solver
14 !> has started to run, allocating additional scratch memory.
15 !===============================================================================
16 
18  USE io_mod
19  USE trace_mod
20  USE patch_mod
21  USE mesh_mod
22  USE link_mod
23  USE task_list_mod
24  USE lagrange_mod
25  USE experiment_mod
26  USE omp_mod
27  USE scaling_mod
28  implicit none
29  private
30  type, public:: rs_t
31  contains
32  procedure, nopass:: init
33  end type
34  character(len=64), save:: meshfile='mhd63.msh'
35  character(len=64), save:: snapfile='mhd63.dat'
36  character(len=64), save:: model='bifrost'
37  integer, save:: order=2
38  integer, save:: verbose=1
39  logical, save:: flip_z = .false.
40  real, save:: b0=0.0
41  type(rs_t), public:: read_snapshot
42 CONTAINS
43 
44 !===============================================================================
45 !===============================================================================
46 SUBROUTINE init (task_list)
47  type(task_list_t), pointer:: task_list
48  integer:: iostat
49  namelist /snapshot_params/ verbose, meshfile, order, verbose, snapfile, model, b0, flip_z
50  !-----------------------------------------------------------------------------
51  rewind(io%input)
52  read (io%input, snapshot_params, iostat=iostat)
53  if (io%master) then
54  write (*, snapshot_params)
55  else
56  verbose = 0
57  end if
58  if (io%restart < 0) then
59  select case(trim(model))
60  case('bifrost')
61  call read_bifrost_snap(task_list)
62  case('stagger')
63  call read_stagger_snap(task_list)
64  case default
65  call io%abort('read_snapshot_mod::init: unknown input model')
66  end select
67  end if
68 END SUBROUTINE init
69 
70 !===============================================================================
71 !> Read a legacy Bifrost raw data file, with auxiliary grid file
72 !===============================================================================
73 SUBROUTINE read_bifrost_snap (task_list)
74  USE kinds_mod
75  type(task_list_t), pointer:: task_list
76  !.............................................................................
77  class(link_t), pointer:: link
78  class(patch_t), pointer:: patch
79  class(mesh_t), pointer, dimension(:):: m
80  integer:: n(3), iv, iiv, recl, ix, iy, iz
81  integer:: zi
82  integer:: mx, my, mz, nx, ny, nz, nv, ii(8)
83  real:: sc_u, sc_b
84  real, dimension(:), allocatable, target:: dxm, dxmdn, xm, xmdn
85  real, dimension(:), allocatable, target:: dym, dymdn, ym, ymdn
86  real, dimension(:), allocatable, target:: dzm, dzmdn, zm, zmdn
87  real, dimension(:,:,:), allocatable:: rbuf, buffer
88  integer, save:: itimer=0
89  integer:: iostat
90  !-----------------------------------------------------------------------------
91  call trace%begin('read_bifrost_snap', itimer=itimer)
92  call io%header('begin ic_mod%read_snapshot: Reading Bifrost legacy data')
93  !-----------------------------------------------------------------------------
94  ! Pick up patch info from task_list
95  !-----------------------------------------------------------------------------
96  patch => task2patch(task_list%head%task)
97  m => patch%mesh
98  nv = 8
99  n = nint(m%b/m%d)
100  !---------------------------------------------------------------------------
101  ! Read the mesh data from the mesh file
102  !---------------------------------------------------------------------------
103  open (io_unit%datain, file=trim(meshfile),status='old', form='formatted')
104  ! X direction
105  read(io_unit%datain, *, iostat=iostat) mx
106  print *, 'mx',mx
107  allocate( dxm(mx), dxmdn(mx), xm(mx+1), xmdn(mx+1))
108  read(io_unit%datain, *, iostat=iostat) (xm(ix),ix=1,mx), (xmdn(ix),ix=1,mx), dxm, dxmdn
109 
110  ! Y direction
111  read(io_unit%datain, *, iostat=iostat) my
112  print *, 'my',my
113  allocate( dym(my), dymdn(my), ym(my+1), ymdn(my+1))
114  read(io_unit%datain, *, iostat=iostat) (ym(iy),iy=1,my), (ymdn(iy),iy=1,my), dym, dymdn
115 
116  ! Z direction
117  read(io_unit%datain, *, iostat=iostat) mz
118  print *, 'mz',mz
119  allocate( dzm(mz), dzmdn(mz), zm(mz ), zmdn(mz ))
120  read(io_unit%datain, *, iostat=iostat) (zm(iz),iz=1,mz), (zmdn(iz),iz=1,mz), dzm, dzmdn
121 
122  if (io%master) then
123  print *,'reading mesh data from ' //trim(meshfile)
124  print *, 'setup dimensions =', n, nv
125  print *, ' file dims =', mx, my, mz
126  print *, ' mesh dims =', m%n, m%gn
127  print *, 'setup dims =', nint(m%b/m%d), m%b
128  end if
129 
130  ! periodic extensions
131  xm(mx+1) = xm(mx) + (xm(mx)-xm(mx-1))
132  ym(my+1) = ym(my) + (ym(my)-ym(my-1))
133  xmdn(mx+1) = xmdn(mx) + (xmdn(mx)-xmdn(mx-1))
134  ymdn(my+1) = ymdn(my) + (ymdn(my)-ymdn(my-1))
135 
136  if (verbose>1) then
137  do iz=1,size(dzm)
138  print *, 'initial_t%stagger: iz, min(dx),min(dy), dz', &
139  iy, minval(dxm), minval(dym), dzm(iz)
140  end do
141  end if
142  if (io%master) then
143  print *, 'min(dx),min(dy),min(dz)', minval(dxm), minval(dym), minval(dzm)
144  print *, mz, zm(1), size(dzm)
145  print *, 'nz for max res:', int((zm(mz)-zm(1))/minval(dzm))+1
146  end if
147  close (io_unit%datain)
148 
149  !---------------------------------------------------------------------------
150  ! Allocate a buffer for reading Bifrost legacy snapshots
151  !---------------------------------------------------------------------------
152  allocate (rbuf(mx,my,mz))
153  if (allocated(rbuf))call io%bits_mem (storage_size(rbuf), mx*my*mz)
154  allocate (buffer(mx+1,my+1,mz))
155  call io%bits_mem (storage_size(buffer), product(shape(buffer)), 'buffer')
156  recl = kind(1.0)*mx*my*mz
157  if (io%master) &
158  print *,'reading MHD data from ',trim(snapfile)
159  open (io_unit%datain, file=trim(snapfile), form='unformatted', access='direct', recl=recl)
160  ! Indices
161  ii=[patch%idx%d,patch%idx%px,patch%idx%py,patch%idx%pz, &
162  patch%idx%e,patch%idx%bx,patch%idx%by,patch%idx%bz]
163  !-----------------------------------------------------------------------------
164  ! Read data and optionally swap z
165  !-----------------------------------------------------------------------------
166  do iv=1,min(nv,8)
167  iiv = ii(iv)
168  if (io%master) print *, 'reading variable index iv =', iv
169  read (io_unit%datain, rec=iv) rbuf(1:mx,1:my,1:mz)
170 
171  print *, 'iv minmax in snapshot', minval(rbuf(1:mx,1:my,1:mz)), maxval(rbuf(1:mx,1:my,1:mz))
172  print *, 'iv minmax location in snapshot', minloc(rbuf(1:mx,1:my,1:mz)), maxloc(rbuf(1:mx,1:my,1:mz))
173 
174  if (flip_z) then
175  do iz=1,mz; do iy=1,my; do ix=1,mx
176  zi = mz -iz + 1
177  buffer(ix,iy,iz) = rbuf(ix,iy,zi)
178  end do; end do; end do
179  else
180  do iz=1,mz; do iy=1,my; do ix=1,mx
181  buffer(ix,iy,iz) = rbuf(ix,iy,iz)
182  end do; end do; end do
183  end if
184  if (verbose == 1) &
185  print *, iv, minval(buffer(1:mx,1:my,1:mz)), &
186  minloc(buffer(1:mx,1:my,1:mz)),'||', &
187  maxval(buffer(1:mx,1:my,1:mz)), &
188  maxloc(buffer(1:mx,1:my,1:mz))
189 
190  !---------------------------------------------------------------------------
191  ! Add periodic extension, if needed
192  !---------------------------------------------------------------------------
193  ! for X
194  if (patch%periodic(1)) then
195  buffer(mx+1,:,:) = buffer(1,:,:)
196  else
197  buffer(mx+1,:,:) = buffer(mx,:,:)
198  end if
199  ! for Y
200  if (patch%periodic(2)) then
201  buffer(:,my+1,:) = buffer(:,1,:)
202  else
203  buffer(:,my+1,:) = buffer(:,my,:)
204  end if
205  ! Z is not expected to be periodic, for now.
206  if (patch%periodic(3)) then
207  call io%abort("read_snapshot_mod::read_bifrost: setup is periodic in Z. That not implemented or tested.")
208  end if
209  !---------------------------------------------------------------------------
210  ! Pick out the piece that belongs to each patch, based on its position
211  !---------------------------------------------------------------------------
212  !$omp parallel default (shared) private(link, patch)
213  !$omp single
214  link => task_list%head
215  do while (associated(link))
216  patch => task2patch(link%task)
217  if (verbose==4) print *, 'read_bifrost_snap spawning task for patch%id =', patch%id
218  !-------------------------------------------------------------------------
219  ! Spawn OMP tasks to interpolate the patch data
220  !-------------------------------------------------------------------------
221  !$omp task default(shared) firstprivate(patch,iiv)
222  call interpolate (iiv, patch, buffer, mx, my, mz, xm, ym, zm, xmdn, ymdn, zmdn)
223  !$omp end task
224  link => link%next
225  end do
226  !$omp end single
227  !$omp end parallel
228  end do
229  close (io_unit%datain)
230  call io%bits_mem (-storage_size(buffer), product(shape(buffer)), 'buffer')
231  deallocate (buffer)
232  if (allocated(rbuf))call io%bits_mem (-storage_size(rbuf), mx*my*mz)
233  if (allocated(rbuf)) deallocate(rbuf)
234  deallocate (xm,ym,zm,xmdn,ymdn,zmdn)
235  !-----------------------------------------------------------------------------
236  ! Reset the vertical magnetic field
237  !-----------------------------------------------------------------------------
238  link => task_list%head
239  do while (associated(link))
240  patch => task2patch(link%task)
241  if (b0 >= 0.0) &
242  patch%mem(:,:,:,patch%idx%bz,1,1) = b0
243  link => link%next
244  end do
245  if (io%master) &
246  print *,'read_bifrost_snap finished'
247  call trace%end(itimer)
248 END SUBROUTINE read_bifrost_snap
249 
250 !===============================================================================
251 !> Read a legacy Stagger raw data file, with auxiliary grid file
252 !===============================================================================
253 SUBROUTINE read_stagger_snap (task_list)
254  type(task_list_t), pointer:: task_list
255  !.............................................................................
256  class(link_t), pointer:: link
257  class(patch_t), pointer:: patch
258  class(mesh_t), pointer, dimension(:):: m
259  integer:: n(3), iv, iiv, recl, ix, iy, iz
260  integer:: zi
261  integer:: mx, my, mz, nx, ny, nz, nv, ii(9), biv
262  real, dimension(:), allocatable, target:: dxm, dxmdn, xm, xmdn
263  real, dimension(:), allocatable, target:: dym, dymdn, ym, ymdn
264  real, dimension(:), allocatable, target:: dzm, dzmdn, zm, zmdn
265  real, dimension(:,:,:), allocatable:: rbuf, buffer
266  integer, save:: itimer=0
267  !-----------------------------------------------------------------------------
268  call trace%begin('read_stagger_snap', itimer=itimer)
269  call io%header('begin ic_mod%read_stagger: Reading legacy Stagger Code data')
270  !-----------------------------------------------------------------------------
271  ! Pick up patch info from task_list
272  !-----------------------------------------------------------------------------
273  patch => task2patch(task_list%head%task)
274  m => patch%mesh
275  nv = patch%nv
276  n = nint(m%b/m%d)
277  n(3) = n(3)+4
278  !---------------------------------------------------------------------------
279  ! Read the mesh meta-data from the .msh file
280  !---------------------------------------------------------------------------
281  open (io_unit%datain, file=trim(meshfile), form='unformatted')
282  read(io_unit%datain) mx, mz, my
283  if (io%master) then
284  print *,'reading mesh data from ',trim(meshfile)
285  print *, 'setup dimensions =', n, nv
286  print *, ' file dims =', mx, mz, my
287  print *, ' mesh dims =', m%n, m%gn
288  print *, 'setup dims =', nint(m%b/m%d), m%b
289  end if
290  ! allocate mesh arrays
291  allocate( dxm(mx), dxmdn(mx), xm(mx+1), xmdn(mx+1))
292  allocate( dym(my), dymdn(my), ym(my+1), ymdn(my+1))
293  allocate( dzm(mz), dzmdn(mz), zm(mz ), zmdn(mz ))
294  ! read meash properties, in order y,z,x -- to make z vertical
295  read(io_unit%datain) dym, dymdn, (ym(iy),iy=1,my), (ymdn(iy),iy=1,my)
296  read(io_unit%datain) dzm, dzmdn, (zm(iz),iz=1,mz), (zmdn(iz),iz=1,mz)
297  read(io_unit%datain) dxm, dxmdn, (xm(ix),ix=1,mx), (xmdn(ix),ix=1,mx)
298  ! periodic extensions
299  xm(mx+1) = xm(mx) + (xm(mx)-xm(mx-1))
300  ym(my+1) = ym(my) + (ym(my)-ym(my-1))
301  xmdn(mx+1) = xmdn(mx) + (xmdn(mx)-xmdn(mx-1))
302  ymdn(my+1) = ymdn(my) + (ymdn(my)-ymdn(my-1))
303  if (verbose>1) then
304  do iz=1,size(dzm)
305  print *, 'initial_t%stagger: iz, min(dx),min(dy), dz', &
306  iy, minval(dxm), minval(dym), dzm(iz)
307  end do
308  end if
309  if (io%master) then
310  print *, 'min(dx),min(dy),min(dz)', minval(dxm), minval(dym), minval(dzm)
311  print *, 'nz for max res:', int((zm(mz)-zm(1))/minval(dzm))+1
312  end if
313  close (io_unit%datain)
314  !---------------------------------------------------------------------------
315  ! Allocate buffers for reading and for transposing to z vertical
316  !---------------------------------------------------------------------------
317  allocate (rbuf(my,mz,mx))
318  call io%bits_mem (storage_size(rbuf), mx*my*mz)
319  allocate (buffer(mx+1,my+1,mz))
320  call io%bits_mem (storage_size(buffer), product(shape(buffer)), 'buffer')
321  recl = 4*product(shape(rbuf)) ! FIXME; this could be compiler dependent
322  if (io%master) &
323  print *,'reading MHD data from ',trim(snapfile)
324  open (io_unit%datain, file=trim(snapfile), form='unformatted', access='direct', recl=recl)
325  !-----------------------------------------------------------------------------
326  ! Read data and swap indices, so z becomes vertical
327  !-----------------------------------------------------------------------------
328  ii = [patch%idx%d, patch%idx%py, patch%idx%pz, patch%idx%px, &
329  patch%idx%e, patch%idx%tt, patch%idx%by, patch%idx%bz, patch%idx%bx]
330  biv=1
331  do iv=1,min(nv,9)
332  if ((iv==6).and.(patch%idx%tt<1)) then
333  biv = biv+1
334  print *, 'skipping temperature', iv
335  !cycle
336  end if
337  iiv = ii(biv)
338  if (io%master) &
339  print *, 'reading variable index iv =', biv,'stored in', ii(biv)
340  read (io_unit%datain, rec=biv) rbuf
341  if (flip_z) then
342  do iz=1,mz; do iy=1,my; do ix=1,mx
343  zi = mz - iz + 1
344  buffer(ix,iy,iz) = rbuf(iy,zi,ix)
345  end do; end do; end do
346  else
347  do iz=1,mz; do iy=1,my; do ix=1,mx
348  buffer(ix,iy,iz) = rbuf(iy,iz,ix)
349  end do; end do; end do
350  end if
351  if (verbose == 1) &
352  print *, biv, minval(rbuf), maxval(rbuf)
353  !---------------------------------------------------------------------------
354  ! Add periodic extension
355  !---------------------------------------------------------------------------
356  buffer(mx+1,:,:) = buffer(1,:,:)
357  buffer(:,my+1,:) = buffer(:,1,:)
358  !---------------------------------------------------------------------------
359  ! Pick out the piece that belongs to each patch, based on its position
360  !---------------------------------------------------------------------------
361  !$omp parallel default (shared) private(link, patch)
362  !$omp single
363  link => task_list%head
364  do while (associated(link))
365  patch => task2patch(link%task)
366  if (verbose==4) print *, 'read_stagger_snap spawning task for patch%id =', patch%id
367  !-------------------------------------------------------------------------
368  ! Spawn OMP tasks to interpolate the patch data
369  !-------------------------------------------------------------------------
370  !$omp task default(shared) firstprivate(patch,iiv)
371  call interpolate (iiv, patch, buffer, mx, my, mz, xm, ym, zm, xmdn, ymdn, zmdn)
372  !$omp end task
373  link => link%next
374  end do
375  !$omp end single
376  !$omp end parallel
377  biv = biv+1
378  end do
379  close (io_unit%datain)
380  deallocate (rbuf, buffer)
381  deallocate (xm,ym,zm,xmdn,ymdn,zmdn)
382  !-----------------------------------------------------------------------------
383  ! Reset the vertical magnetic field
384  !-----------------------------------------------------------------------------
385  link => task_list%head
386  do while (associated(link))
387  patch => task2patch(link%task)
388  if (b0 >= 0.0) &
389  patch%mem(:,:,:,patch%idx%bz,1,1) = b0
390  link => link%next
391  end do
392  if (io%master) &
393  print *,'read_stagger_snap finished'
394  call trace%end(itimer)
395 END SUBROUTINE read_stagger_snap
396 
397 !===============================================================================
398 !> Thread safe interpolation work, where only iv, patch, ff changes in the calls
399 !===============================================================================
400 SUBROUTINE interpolate (iv, patch, buffer, mx, my, mz, xm, ym, zm, xmdn, ymdn, zmdn)
401  class(patch_t), pointer:: patch
402  real, dimension(:,:,:), allocatable:: buffer
403  integer:: iv, mx, my, mz
404  real, dimension(:), allocatable, target:: dxm, dxmdn, xm, xmdn
405  real, dimension(:), allocatable, target:: dym, dymdn, ym, ymdn
406  real, dimension(:), allocatable, target:: dzm, dzmdn, zm, zmdn
407  !.............................................................................
408  class(mesh_t), pointer, dimension(:):: m
409  integer:: nx, ny, nz, o(3)
410  real, dimension(:), pointer:: xi, yi, zi, xo, yo, zo
411  !-----------------------------------------------------------------------------
412  m => patch%mesh
413  nx = m(1)%gn
414  ny = m(2)%gn
415  nz = m(3)%gn
416  allocate (xo(nx), yo(ny), zo(nz))
417  associate(m1=>m(1), m2=>m(2), m3=>m(3))
418  xo = m1%p + m1%r + m1%h(iv)*m1%d
419  yo = m2%p + m2%r + m2%h(iv)*m2%d
420  zo = m3%p + m3%r + m3%h(iv)*m3%d
421  end associate
422  xi => xm
423  yi => ym
424  zi => zm
425  if (iv==patch%idx%px.or.iv==patch%idx%bx) xi => xmdn
426  if (iv==patch%idx%py.or.iv==patch%idx%by) yi => ymdn
427  if (iv==patch%idx%pz.or.iv==patch%idx%bz) zi => zmdn
428  xo = modulo(xo-xi(1),xm(mx+1)-xm(1))
429  yo = modulo(yo-yi(1),ym(my+1)-ym(1))
430  associate(ff => patch%mem(:,:,:,iv,patch%it,1))
431  if (order==1) then
432  call lagrange%trilinear3d (xi, yi, zi, buffer, &
433  xo, yo, zo, ff, order)
434  else
435  call lagrange%interpolate3d (xi, yi, zi, buffer, &
436  xo, yo, zo, ff, order)
437  end if
438  !-----------------------------------------------------------------------
439  ! Diagnostic output
440  !-----------------------------------------------------------------------
441  if (verbose==2) then
442  print '(i5,i3,2x,2(2x,a,3(2x,2f7.2)))', patch%id, iv, &
443  'xyz(in):', &
444  minval(xi), maxval(xi), minval(yi), maxval(yi), minval(zi), maxval(zi), &
445  'xyz(out):', &
446  minval(xo), maxval(xo), minval(yo), maxval(yo), minval(zo), maxval(zo)
447  else if (verbose==3) then
448  o = maxloc(ff)
449  print '(1x,a,2i3,3i4,2x,3f8.3,1p,2(2x,2e12.4))', &
450  'omp,iv,ix,iy,iz,pos,bufmin,bufmax,fmin,fmax =', omp%thread, iv, o, m%p, &
451  minval(buffer), maxval(buffer), &
452  minval( ff), maxval( ff)
453  end if
454  end associate
455  deallocate (xo, yo, zo)
456 END SUBROUTINE
457 
458 END MODULE read_snapshot_mod
Define code units, in terms of (here) CGS units.
Definition: scaling_mod.f90:4
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
Module for Lagrange interpolation.
Definition: lagrange_mod.f90:4
Task list data type, with methods for startup and updates. Message handling is inherited from the tas...
Template module for mesh.
Definition: mesh_mod.f90:4
Definition: io_mod.f90:4
IC module for reading in legacy Stagger snapshots into existing patches.