32 procedure,
nopass:: init
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.
41 type(
rs_t),
public:: read_snapshot
46 SUBROUTINE init (task_list)
49 namelist /snapshot_params/ verbose, meshfile, order, verbose, snapfile, model, b0, flip_z
52 read (io%input, snapshot_params, iostat=iostat)
54 write (*, snapshot_params)
58 if (io%restart < 0)
then 59 select case(trim(model))
61 call read_bifrost_snap(task_list)
63 call read_stagger_snap(task_list)
65 call io%abort(
'read_snapshot_mod::init: unknown input model')
73 SUBROUTINE read_bifrost_snap (task_list)
77 class(
link_t),
pointer:: link
79 class(
mesh_t),
pointer,
dimension(:):: m
80 integer:: n(3), iv, iiv, recl, ix, iy, iz
82 integer:: mx, my, mz, nx, ny, nz, nv, ii(8)
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
91 call trace%begin(
'read_bifrost_snap', itimer=itimer)
92 call io%header(
'begin ic_mod%read_snapshot: Reading Bifrost legacy data')
96 patch => task2patch(task_list%head%task)
103 open (io_unit%datain, file=trim(meshfile),status=
'old', form=
'formatted')
105 read(io_unit%datain, *, iostat=iostat) 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
111 read(io_unit%datain, *, iostat=iostat) 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
117 read(io_unit%datain, *, iostat=iostat) 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
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
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))
138 print *,
'initial_t%stagger: iz, min(dx),min(dy), dz', &
139 iy, minval(dxm), minval(dym), dzm(iz)
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
147 close (io_unit%datain)
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
158 print *,
'reading MHD data from ',trim(snapfile)
159 open (io_unit%datain, file=trim(snapfile), form=
'unformatted', access=
'direct', recl=recl)
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]
168 if (io%master) print *,
'reading variable index iv =', iv
169 read (io_unit%datain, rec=iv) rbuf(1:mx,1:my,1:mz)
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))
175 do iz=1,mz;
do iy=1,my;
do ix=1,mx
177 buffer(ix,iy,iz) = rbuf(ix,iy,zi)
178 end do; end do; end do
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
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))
194 if (patch%periodic(1))
then 195 buffer(mx+1,:,:) = buffer(1,:,:)
197 buffer(mx+1,:,:) = buffer(mx,:,:)
200 if (patch%periodic(2))
then 201 buffer(:,my+1,:) = buffer(:,1,:)
203 buffer(:,my+1,:) = buffer(:,my,:)
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.")
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
222 call interpolate (iiv, patch, buffer, mx, my, mz, xm, ym, zm, xmdn, ymdn, zmdn)
229 close (io_unit%datain)
230 call io%bits_mem (-storage_size(buffer), product(shape(buffer)),
'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)
238 link => task_list%head
239 do while (
associated(link))
240 patch => task2patch(link%task)
242 patch%mem(:,:,:,patch%idx%bz,1,1) = b0
246 print *,
'read_bifrost_snap finished' 247 call trace%end(itimer)
248 END SUBROUTINE read_bifrost_snap
253 SUBROUTINE read_stagger_snap (task_list)
254 type(task_list_t),
pointer:: task_list
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
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
268 call trace%begin(
'read_stagger_snap', itimer=itimer)
269 call io%header(
'begin ic_mod%read_stagger: Reading legacy Stagger Code data')
273 patch => task2patch(task_list%head%task)
281 open (io_unit%datain, file=trim(meshfile), form=
'unformatted')
282 read(io_unit%datain) mx, mz, my
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
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 ))
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)
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))
305 print *,
'initial_t%stagger: iz, min(dx),min(dy), dz', &
306 iy, minval(dxm), minval(dym), dzm(iz)
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
313 close (io_unit%datain)
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))
323 print *,
'reading MHD data from ',trim(snapfile)
324 open (io_unit%datain, file=trim(snapfile), form=
'unformatted', access=
'direct', recl=recl)
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]
332 if ((iv==6).and.(patch%idx%tt<1))
then 334 print *,
'skipping temperature', iv
339 print *,
'reading variable index iv =', biv,
'stored in', ii(biv)
340 read (io_unit%datain, rec=biv) rbuf
342 do iz=1,mz;
do iy=1,my;
do ix=1,mx
344 buffer(ix,iy,iz) = rbuf(iy,zi,ix)
345 end do; end do; end do
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
352 print *, biv, minval(rbuf), maxval(rbuf)
356 buffer(mx+1,:,:) = buffer(1,:,:)
357 buffer(:,my+1,:) = buffer(:,1,:)
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
371 call interpolate (iiv, patch, buffer, mx, my, mz, xm, ym, zm, xmdn, ymdn, zmdn)
379 close (io_unit%datain)
380 deallocate (rbuf, buffer)
381 deallocate (xm,ym,zm,xmdn,ymdn,zmdn)
385 link => task_list%head
386 do while (
associated(link))
387 patch => task2patch(link%task)
389 patch%mem(:,:,:,patch%idx%bz,1,1) = b0
393 print *,
'read_stagger_snap finished' 394 call trace%end(itimer)
395 END SUBROUTINE read_stagger_snap
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
408 class(mesh_t),
pointer,
dimension(:):: m
409 integer:: nx, ny, nz, o(3)
410 real,
dimension(:),
pointer:: xi, yi, zi, xo, yo, zo
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
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))
432 call lagrange%trilinear3d (xi, yi, zi, buffer, &
433 xo, yo, zo, ff, order)
435 call lagrange%interpolate3d (xi, yi, zi, buffer, &
436 xo, yo, zo, ff, order)
442 print
'(i5,i3,2x,2(2x,a,3(2x,2f7.2)))', patch%id, iv, &
444 minval(xi), maxval(xi), minval(yi), maxval(yi), minval(zi), maxval(zi), &
446 minval(xo), maxval(xo), minval(yo), maxval(yo), minval(zo), maxval(zo)
447 else if (verbose==3)
then 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)
455 deallocate (xo, yo, zo)
Define code units, in terms of (here) CGS units.
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Module for Lagrange interpolation.
Task list data type, with methods for startup and updates. Message handling is inherited from the tas...
Template module for mesh.
Module with list handling for generic class task_t objects.
IC module for reading in legacy Stagger snapshots into existing patches.