MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
convert.t
Go to the documentation of this file.
1 !=============================================================================
2 subroutine generate_plotfile
8 !-----------------------------------------------------------------------------
9 
10 if(mype==0.and.level_io>0) write(unitterm,*)'reset tree to fixed level=',level_io
11 if(level_io>0 .or. level_io_min.ne.1 .or. level_io_max.ne.nlevelshi) then
13 else if(.not. phys_req_diagonal) then
14  call getbc(global_time,0.d0,ps,1,nwflux+nwaux)
15 end if
16 
17 select case(convert_type)
18  case('tecplot','tecplotCC','tecline')
19  call tecplot(unitconvert)
20  case('tecplotmpi','tecplotCCmpi','teclinempi')
22  case('vtu','vtuCC')
24  case('vtumpi','vtuCCmpi')
26  case('vtuB','vtuBCC','vtuBmpi','vtuBCCmpi')
28  case('pvtumpi','pvtuCCmpi')
30  case('pvtuBmpi','pvtuBCCmpi')
32  case('vtimpi','vtiCCmpi')
34  case('onegrid','onegridmpi')
35  call onegrid(unitconvert)
36  case('oneblock','oneblockB')
37  call oneblock(unitconvert)
38  case('user','usermpi')
39  if (.not. associated(usr_special_convert)) then
40  call mpistop("usr_special_convert not defined")
41  else
43  end if
44  case default
45  call mpistop("Error in generate_plotfile: Unknown convert_type")
46 end select
47 
48 
49 ! output synthetic euv emission
50 if (ndim==3) then
53 endif
54 
55 end subroutine generate_plotfile
56 !=============================================================================
57 subroutine oneblock(qunit)
58 ! this is for turning an AMR run into a single block
59 ! the data will be all on selected level level_io
60 
61 ! this version should work for any dimension
62 ! only writes w_write selected 1:nw variables, also nwauxio
63 ! may use saveprim to switch to primitives
64 ! this version can not work on multiple CPUs
65 ! does not renormalize variables
66 
67 ! header info differs from onegrid below
68 
69 ! ASCII or binary output
70 
74 use mod_physics
76 integer, intent(in) :: qunit
77 
78 integer :: Morton_no,igrid,ix^D,ig^D,level
79 integer, pointer :: ig_to_igrid(:^d&,:)
80 logical :: fileopen,writeblk(max_blocks)
81 character(len=80) :: filename
82 integer :: filenr,ncells,ncells^D,ncellg,ncellx^D,jg^D,jig^D
83 
84 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
85 character(len=1024) :: outfilehead
86 
87 double precision :: wval1,xval1
88 double precision, dimension({^D&1:1},1:nw+nwauxio) :: wval
89 double precision, dimension({^D&1:1},1:ndim) :: xval
90 double precision:: normconv(0:nw+nwauxio)
91 
92 integer :: iw,iiw,writenw,iwrite(1:nw+nwauxio),iigrid,idim
93 logical :: patchw(ixg^t)
94 !-----------------------------------------------------------------------------
95 
96 if(level_io<1)then
97  call mpistop('please specify level_io>0 for usage with oneblock')
98 end if
99 
100 if(npe>1)then
101  if(mype==0) print *,'ONEBLOCK as yet to be parallelized'
102  call mpistop('npe>1, oneblock')
103 end if
104 
105 ! only variables selected by w_write will be written out
106 normconv(0:nw+nwauxio)=one
107 normconv(0) = length_convert_factor
108 normconv(1:nw) = w_convert_factor
109 writenw=count(w_write(1:nw))+nwauxio
110 iiw=0
111 do iw =1,nw
112  if (.not.w_write(iw))cycle
113  iiw=iiw+1
114  iwrite(iiw)=iw
115 end do
116 if(nwauxio>0)then
117  do iw =nw+1,nw+nwauxio
118  iiw=iiw+1
119  iwrite(iiw)=iw
120  end do
121 endif
122 
123 allocate(ig_to_igrid(ng^d(level_io),0:npe-1))
124 ig_to_igrid=-1
125 writeblk=.false.
126 do morton_no=morton_start(mype),morton_stop(mype)
127  igrid=sfc_to_igrid(morton_no)
128  level=node(plevel_,igrid)
129  ig^d=igrid_to_node(igrid,mype)%node%ig^d;
130  ig_to_igrid(ig^d,mype)=igrid
131  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
132  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
133  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
134  writeblk(igrid)=.true.
135  end if
136 end do
137 
138 call getheadernames(wnamei,xandwnamei,outfilehead)
139 ncells=0
140 ncells^d=0;
141 ncellg=(^d&(ixmhi^d-ixmlo^d+1)*)
142 ncellx^d=ixmhi^d-ixmlo^d+1\
143 {do ig^d=1,ng^d(level_io)\}
144  igrid=ig_to_igrid(ig^d,mype)
145  if(writeblk(igrid)) go to 20
146 {end do\}
147 20 continue
148 jg^d=ig^d;
149 {
150 jig^dd=jg^dd;
151 do ig^d=1,ng^d(level_io)
152  jig^d=ig^d
153  igrid=ig_to_igrid(jig^dd,mype)
154  if(writeblk(igrid)) ncells^d=ncells^d+ncellx^d
155 end do
156 \}
157 
158 do iigrid=1,igridstail; igrid=igrids(iigrid)
159  if(.not.writeblk(igrid)) cycle
160  ncells=ncells+ncellg
161  ps1(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
162 
163  if (nwauxio > 0) then
164  if (.not. associated(usr_aux_output)) then
165  call mpistop("usr_aux_output not defined")
166  else
167  call usr_aux_output(ixg^ll,ixm^ll^ladd1, &
168  ps1(igrid)%w,ps(igrid)%x,normconv)
169  end if
170  end if
171 end do
172 
173 if (saveprim) then
174  do iigrid=1,igridstail; igrid=igrids(iigrid)
175  if (.not.writeblk(igrid)) cycle
176  call phys_to_primitive(ixg^ll,ixg^ll^lsub1,ps1(igrid)%w,ps(igrid)%x)
177  if(b0field) then
178  ! add background magnetic field B0 to B
179  ps1(igrid)%w(ixg^t,iw_mag(:))=ps1(igrid)%w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
180  end if
181  end do
182 else
183  do iigrid=1,igridstail; igrid=igrids(iigrid)
184  if (.not.writeblk(igrid)) cycle
185  if (b0field) then
186  ! add background magnetic field B0 to B
187  if(phys_energy) &
188  ps1(igrid)%w(ixg^t,iw_e)=ps1(igrid)%w(ixg^t,iw_e)+0.5d0*sum(ps(igrid)%B0(ixg^t,:,0)**2,dim=ndim+1) &
189  + sum(ps1(igrid)%w(ixg^t,iw_mag(:))*ps(igrid)%B0(ixg^t,:,0),dim=ndim+1)
190  ps1(igrid)%w(ixg^t,iw_mag(:))=ps1(igrid)%w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
191  end if
192  end do
193 end if
194 
195 master_cpu_open : if (mype == 0) then
196  inquire(qunit,opened=fileopen)
197  if (.not.fileopen) then
198  ! generate filename
199  filenr=snapshotini
200  if (autoconvert) filenr=snapshotnext
201  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".blk"
202  select case(convert_type)
203  case("oneblock")
204  open(qunit,file=filename,status='unknown')
205  write(qunit,*) trim(outfilehead)
206  write(qunit,*) ncells,ncells^d
207  write(qunit,*) global_time*time_convert_factor
208  case("oneblockB")
209  open(qunit,file=filename,form='unformatted',status='unknown')
210  write(qunit) outfilehead
211  write(qunit) ncells,ncells^d
212  write(qunit) global_time*time_convert_factor
213  end select
214  end if
215 end if master_cpu_open
216 
217 {^ifthreed
218 do ig3=1,ng3(level_io)
219  do ix3=ixmlo3,ixmhi3}
220 
221  {^nooned
222  do ig2=1,ng2(level_io)
223  do ix2=ixmlo2,ixmhi2}
224 
225  do ig1=1,ng1(level_io)
226  igrid=ig_to_igrid(ig^d,mype)
227  if(.not.writeblk(igrid)) cycle
228  do ix1=ixmlo1,ixmhi1
229  master_write : if(mype==0) then
230  select case(convert_type)
231  case("oneblock")
232  write(qunit,fmt="(100(e14.6))") &
233  ps(igrid)%x(ix^d,1:ndim)*normconv(0),&
234  (ps1(igrid)%w(ix^d,iwrite(iw))*normconv(iwrite(iw)),iw=1,writenw)
235  case("oneblockB")
236  write(qunit) real(ps(igrid)%x(ix^D,1:ndim)*normconv(0)),&
237  (real(ps1(igrid)%w(ix^D,iwrite(iw))*normconv(iwrite(iw))),iw=1,writenw)
238  end select
239  end if master_write
240  end do
241  end do
242  {^nooned
243  end do
244  end do}
245  {^ifthreed
246  end do
247 end do}
248 
249 close(qunit)
250 
251 end subroutine oneblock
252 !=============================================================================
253 subroutine onegrid(qunit)
255 ! this is for turning an AMR run into a single grid
256 ! this version should work for any dimension, can be in parallel
257 ! in 1D, should behave much like oneblock, except for header info
258 
259 ! only writes all 1:nw variables, no nwauxio
260 ! may use saveprim to switch to primitives
261 ! this version can work on multiple CPUs
262 ! does not renormalize variables
263 ! ASCII output
264 
267 use mod_physics
269 integer, intent(in) :: qunit
270 
271 integer :: Morton_no,igrid,ix^D,iw
272 logical :: fileopen
273 character(len=80) :: filename
274 integer :: filenr
275 
276 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
277 character(len=1024) :: outfilehead
278 
279 !.. MPI variables ..
280 integer :: igrid_recv,ipe
281 double precision :: w_recv(ixg^t,1:nw),x_recv(ixg^t,1:ndim)
282 integer, allocatable :: intstatus(:,:)
283 
284 !-----------------------------------------------------------------------------
285 
286 if(nwauxio>0)then
287  if(mype==0) print *,'ONEGRID to be used without nwauxio'
288  call mpistop('nwauxio>0, onegrid')
289 end if
290 
291 if(saveprim)then
292  if(mype==0.and.nwaux>0) print *,'warning: ONEGRID used with saveprim, check auxiliaries'
293 end if
294 
295 
296 
297 master_cpu_open : if (mype == 0) then
298  call getheadernames(wnamei,xandwnamei,outfilehead)
299  write(outfilehead,'(a)') "#"//" "//trim(outfilehead)
300  inquire(qunit,opened=fileopen)
301  if (.not.fileopen) then
302  ! generate filename
303  filenr=snapshotini
304  if (autoconvert) filenr=snapshotnext
305  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".blk"
306  open(qunit,file=filename,status='unknown')
307  end if
308  write(qunit,"(a)")outfilehead
309  write(qunit,"(i7)") ( {^d&(ixmhi^d-ixmlo^d+1)*} )*(morton_stop(npe-1)-morton_start(0)+1)
310 end if master_cpu_open
311 
312 do morton_no=morton_start(mype),morton_stop(mype)
313  igrid=sfc_to_igrid(morton_no)
314  if(saveprim) call phys_to_primitive(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
315  if (mype/=0)then
316  itag=morton_no
317  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
318  call mpi_send(ps(igrid)%x,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
319  itag=igrid
320  call mpi_send(ps(igrid)%w,1,type_block_io, 0,itag,icomm,ierrmpi)
321  else
322  {do ix^db=ixmlo^db,ixmhi^db\}
323  do iw=1,nw
324  if( dabs(ps(igrid)%w(ix^d,iw)) < 1.0d-32 ) ps(igrid)%w(ix^d,iw) = zero
325  enddo
326  write(qunit,fmt="(100(e14.6))") ps(igrid)%x(ix^d,1:ndim)&
327  ,ps(igrid)%w(ix^d,1:nw)
328  {end do\}
329  end if
330 end do
331 
332 if(mype==0.and.npe>1) allocate(intstatus(mpi_status_size,1))
333 
334 manycpu : if (npe>1) then
335  if (mype==0) then
336  loop_cpu : do ipe =1, npe-1
337  loop_morton : do morton_no=morton_start(ipe),morton_stop(ipe)
338  itag=morton_no
339  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
340  call mpi_recv(x_recv,1,type_block_xcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
341  itag=igrid_recv
342  call mpi_recv(w_recv,1,type_block_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
343  {do ix^db=ixmlo^db,ixmhi^db\}
344  do iw=1,nw
345  if( dabs(ps(igrid)%w(ix^d,iw)) < smalldouble ) ps(igrid)%w(ix^d,iw) = zero
346  enddo
347  write(qunit,fmt="(100(e14.6))") x_recv(ix^d,1:ndim)&
348  ,w_recv(ix^d,1:nw)
349  {end do\}
350  end do loop_morton
351  end do loop_cpu
352  end if
353 end if manycpu
354 
355 if (npe>1) then
356  call mpi_barrier(icomm,ierrmpi)
357  if(mype==0)deallocate(intstatus)
358 endif
359 
360 if(mype==0) close(qunit)
361 end subroutine onegrid
362 !============================================================================
363 subroutine tecplot(qunit)
365 ! output for tecplot (ASCII format)
366 ! not parallel, uses calc_grid to compute nwauxio variables
367 ! allows renormalizing using convert factors
368 
371 
372 integer, intent(in) :: qunit
373 
374 integer:: igrid,iigrid,level,igonlevel,iw,idim,ix^D
375 integer:: NumGridsOnLevel(1:nlevelshi)
376 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,ixC^L,ixCC^L
377 
378 integer :: nodes, elems
379 
380 
381 double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
382 
383 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
384 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
385 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
386 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
387 
388 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
389 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
390 double precision, dimension(0:nw+nwauxio) :: normconv
391 logical :: fileopen,first
392 character(len=80) :: filename
393 integer :: filenr
394 
395 !!! possible length conflict
396 character(len=1024) :: tecplothead
397 
398 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
399 character(len=1024) :: outfilehead
400 !-----------------------------------------------------------------------------
401 
402 if(npe>1)then
403  if(mype==0) print *,'tecplot not parallel, use tecplotmpi'
404  call mpistop('npe>1, tecplot')
405 end if
406 
407 if(nw/=count(w_write(1:nw)))then
408  if(mype==0) print *,'tecplot does not use w_write=F'
409  call mpistop('w_write, tecplot')
410 end if
411 
412 if(nocartesian)then
413  if(mype==0) print *,'tecplot with nocartesian'
414 endif
415 
416 inquire(qunit,opened=fileopen)
417 if (.not.fileopen) then
418  ! generate filename
419  filenr=snapshotini
420  if (autoconvert) filenr=snapshotnext
421  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".plt"
422  open(qunit,file=filename,status='unknown')
423 end if
424 
425 call getheadernames(wnamei,xandwnamei,outfilehead)
426 
427 write(tecplothead,'(a)') "VARIABLES = "//trim(outfilehead)
428 
429 write(qunit,'(a)') tecplothead(1:len_trim(tecplothead))
430 
431 numgridsonlevel(1:nlevelshi)=0
432 do level=levmin,levmax
433  numgridsonlevel(level)=0
434  do iigrid=1,igridstail; igrid=igrids(iigrid);
435  if (node(plevel_,igrid)/=level) cycle
436  numgridsonlevel(level)=numgridsonlevel(level)+1
437  end do
438 end do
439 
440 nx^d=ixmhi^d-ixmlo^d+1;
441 nxc^d=nx^d+1;
442 
443 {^ifoned
444 if(convert_type=='tecline') then
445  nodes=0
446  elems=0
447  do level=levmin,levmax
448  nodes=nodes + numgridsonlevel(level)*{nxc^d*}
449  elems=elems + numgridsonlevel(level)*{nx^d*}
450  enddo
451 
452  write(qunit,"(a,i7,a,1pe12.5,a)") &
453  'ZONE T="all levels", I=',elems, &
454  ', SOLUTIONTIME=',global_time*time_convert_factor,', F=POINT'
455 
456  igonlevel=0
457  do iigrid=1,igridstail; igrid=igrids(iigrid);
458  call calc_x(igrid,xc,xcc)
459  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,ixc^l,ixcc^l,.true.)
460  {do ix^db=ixccmin^db,ixccmax^db\}
461  x_tec(1:ndim)=xcc_tmp(ix^d,1:ndim)*normconv(0)
462  w_tec(1:nw+nwauxio)=wcc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
463  !write(qunit,fmt="(100(e14.6))") x_TEC, w_TEC
464  write(qunit,fmt="(100(e24.16))") x_tec, w_tec
465  {end do\}
466  enddo
467  close(qunit)
468 else
469 }
470 do level=levmin,levmax
471  nodesonlevel=numgridsonlevel(level)*{nxc^d*}
472  elemsonlevel=numgridsonlevel(level)*{nx^d*}
473  ! for all tecplot variants coded up here, we let the TECPLOT ZONES coincide
474  ! with the AMR grid LEVEL. Other options would be
475  ! let each grid define a zone: inefficient for TECPLOT internal workings
476  ! hence not implemented
477  ! let entire octree define 1 zone: no difference in interpolation
478  ! properties across TECPLOT zones detected as yet, hence not done
479  select case(convert_type)
480  case('tecplot')
481  ! in this option, we store the corner coordinates, as well as the corner
482  ! values of all variables (obtained by averaging). This allows POINT packaging,
483  ! and thus we can save full grid info by using one call to calc_grid
484  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
485  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
486  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
487  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
488  do iigrid=1,igridstail; igrid=igrids(iigrid);
489  if (node(plevel_,igrid)/=level) cycle
490  call calc_x(igrid,xc,xcc)
491  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
492  ixc^l,ixcc^l,.true.)
493  {do ix^db=ixcmin^db,ixcmax^db\}
494  x_tec(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0)
495  w_tec(1:nw+nwauxio)=wc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
496  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
497  {end do\}
498  enddo
499  case('tecplotCC')
500  ! in this option, we store the corner coordinates, and the cell center
501  ! values of all variables. Due to this mix of corner/cell center, we must
502  ! use BLOCK packaging, and thus we have enormous overhead by using
503  ! calc_grid repeatedly to merely fill values of cell corner coordinates
504  ! and cell center values per dimension, per variable
505  if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
506  if(nw+nwauxio==1)then
507  ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
508  ! and just set [ndim+1]
509  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
510  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
511  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
512  ndim+1,']=CELLCENTERED), ZONETYPE=', &
513  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
514  else
515  if(ndim+nw+nwauxio<10) then
516  ! difference only in length of integer format specification for ndim+nw+nwauxio
517  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
518  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
519  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
520  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
521  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
522  else
523  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
524  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
525  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
526  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
527  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
528  endif
529  endif
530  do idim=1,ndim
531  first=(idim==1)
532  do iigrid=1,igridstail; igrid=igrids(iigrid);
533  if (node(plevel_,igrid)/=level) cycle
534  call calc_x(igrid,xc,xcc)
535  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
536  ixc^l,ixcc^l,first)
537  write(qunit,fmt="(100(e14.6))") xc_tmp(ixc^s,idim)*normconv(0)
538  enddo
539  enddo
540  do iw=1,nw+nwauxio
541  do iigrid=1,igridstail; igrid=igrids(iigrid);
542  if (node(plevel_,igrid)/=level) cycle
543  call calc_x(igrid,xc,xcc)
544  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
545  ixc^l,ixcc^l,.true.)
546  write(qunit,fmt="(100(e14.6))") wcc_tmp(ixcc^s,iw)*normconv(iw)
547  enddo
548  enddo
549  case default
550  call mpistop('no such tecplot type')
551  end select
552  igonlevel=0
553  do iigrid=1,igridstail; igrid=igrids(iigrid);
554  if (node(plevel_,igrid)/=level) cycle
555  igonlevel=igonlevel+1
556  call save_conntec(qunit,igrid,igonlevel)
557  end do
558 end do
559 {^ifoned endif}
560 
561 close(qunit)
562 
563 end subroutine tecplot
564 !=============================================================================
565 subroutine save_conntec(qunit,igrid,igonlevel)
567 ! this saves the basic line, quad and brick connectivity,
568 ! as used by TECPLOT file outputs for unstructured grid
569 
571 
572 integer, intent(in) :: qunit, igrid, igonlevel
573 
574 integer :: nx^D, nxC^D, ix^D
575 {^ifoned integer, external:: nodenumbertec1D \}
576 {^iftwod integer, external:: nodenumbertec2D \}
577 {^ifthreed integer, external:: nodenumbertec3D \}
578 !-----------------------------------------------------------------------------
579 nx^d=ixmhi^d-ixmlo^d+1;
580 nxc^d=nx^d+1;
581 
582 ! connectivity list
583 {do ix^db=1,nx^db\}
584  {^ifthreed
585  ! basic brick connectivity
586  write(qunit,'(8(i7,1x))') &
587  nodenumbertec3d(ix1, ix2-1,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
588  nodenumbertec3d(ix1+1,ix2-1,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
589  nodenumbertec3d(ix1+1,ix2 ,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
590  nodenumbertec3d(ix1 ,ix2 ,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
591  nodenumbertec3d(ix1 ,ix2-1,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
592  nodenumbertec3d(ix1+1,ix2-1,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
593  nodenumbertec3d(ix1+1,ix2 ,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
594  nodenumbertec3d(ix1 ,ix2 ,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid)
595  }
596  {^iftwod
597  ! basic quadrilateral connectivity
598  write(qunit,'(4(i7,1x))') &
599  nodenumbertec2d(ix1, ix2-1,nxc1,nxc2,igonlevel,igrid),&
600  nodenumbertec2d(ix1+1,ix2-1,nxc1,nxc2,igonlevel,igrid),&
601  nodenumbertec2d(ix1+1,ix2 ,nxc1,nxc2,igonlevel,igrid),&
602  nodenumbertec2d(ix1 ,ix2 ,nxc1,nxc2,igonlevel,igrid)
603  }
604  {^ifoned
605  ! basic line connectivity
606  write(qunit,'(2(i7,1x))') nodenumbertec1d(ix1,nxc1,igonlevel,igrid),&
607  nodenumbertec1d(ix1+1,nxc1,igonlevel,igrid)
608  }
609 {end do\}
610 
611 end subroutine save_conntec
612 !=============================================================================
613 integer function nodenumbertec1d(i1,nx1,ig,igrid)
615 integer, intent(in):: i1,nx1,ig,igrid
616 !-----------------------------------------------------------------------------
617 nodenumbertec1d=i1+(ig-1)*nx1
618 
619 if(nodenumbertec1d>9999999)call mpistop("too large nodenumber")
620 end function nodenumbertec1d
621 !====================================================================================
622 integer function nodenumbertec2d(i1,i2,nx1,nx2,ig,igrid)
624 integer, intent(in):: i1,i2,nx1,nx2,ig,igrid
625 !-----------------------------------------------------------------------------
626 nodenumbertec2d=i1+i2*nx1+(ig-1)*nx1*nx2
627 
628 if(nodenumbertec2d>9999999)call mpistop("too large nodenumber")
629 end function nodenumbertec2d
630 !====================================================================================
631 integer function nodenumbertec3d(i1,i2,i3,nx1,nx2,nx3,ig,igrid)
633 integer, intent(in):: i1,i2,i3,nx1,nx2,nx3,ig,igrid
634 !-----------------------------------------------------------------------------
635 nodenumbertec3d=i1+i2*nx1+i3*nx1*nx2+(ig-1)*nx1*nx2*nx3
636 
637 if(nodenumbertec3d>9999999)call mpistop("too large nodenumber")
638 end function nodenumbertec3d
639 !=============================================================================
640 subroutine unstructuredvtk(qunit)
642 ! output for vtu format to paraview
643 ! not parallel, uses calc_grid to compute nwauxio variables
644 ! allows renormalizing using convert factors
645 
648 
649 integer, intent(in) :: qunit
650 
651 double precision :: x_VTK(1:3)
652 
653 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
654 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
655 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
656 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
657 
658 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
659 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
660 double precision, dimension(0:nw+nwauxio) :: normconv
661 integer:: igrid,iigrid,level,igonlevel,icel,ixC^L,ixCC^L,iw
662 integer:: NumGridsOnLevel(1:nlevelshi)
663 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,VTK_type,ix^D
664 
665 character(len=80):: filename
666 integer :: filenr
667 
668 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
669 character(len=1024) :: outfilehead
670 
671 logical :: fileopen
672 !-----------------------------------------------------------------------------
673 
674 if(npe>1)then
675  if(mype==0) print *,'unstructuredvtk not parallel, use vtumpi'
676  call mpistop('npe>1, unstructuredvtk')
677 end if
678 
679 inquire(qunit,opened=fileopen)
680 if(.not.fileopen)then
681  ! generate filename
682  filenr=snapshotini
683  if (autoconvert) filenr=snapshotnext
684  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
685  ! Open the file for the header part
686  open(qunit,file=filename,status='unknown')
687 endif
688 
689 call getheadernames(wnamei,xandwnamei,outfilehead)
690 
691 ! generate xml header
692 write(qunit,'(a)')'<?xml version="1.0"?>'
693 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
694 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
695 write(qunit,'(a)')'<UnstructuredGrid>'
696 write(qunit,'(a)')'<FieldData>'
697 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
698  'NumberOfTuples="1" format="ascii">'
699 write(qunit,*) dble(dble(global_time)*time_convert_factor)
700 write(qunit,'(a)')'</DataArray>'
701 write(qunit,'(a)')'</FieldData>'
702 
703 ! number of cells, number of corner points, per grid.
704 nx^d=ixmhi^d-ixmlo^d+1;
705 nxc^d=nx^d+1;
706 nc={nx^d*}
707 np={nxc^d*}
708 
709 ! Note: using the w_write, writelevel, writespshift
710 ! we can clip parts of the grid away, select variables, levels etc.
711 do level=levmin,levmax
712  if (writelevel(level)) then
713  do iigrid=1,igridstail; igrid=igrids(iigrid);
714  if (node(plevel_,igrid)/=level) cycle
715  ! only output a grid when fully within clipped region selected
716  ! by writespshift array
717  if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
718  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
719  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
720  call calc_x(igrid,xc,xcc)
721  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
722  ixc^l,ixcc^l,.true.)
723  select case(convert_type)
724  case('vtu')
725  ! we write out every grid as one VTK PIECE
726  write(qunit,'(a,i7,a,i7,a)') &
727  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
728  write(qunit,'(a)')'<PointData>'
729  do iw=1,nw+nwauxio
730  if(iw<=nw) then
731  if(.not.w_write(iw)) cycle
732  endif
733 
734  write(qunit,'(a,a,a)')&
735  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
736  write(qunit,'(200(1pe14.6))') {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
737  write(qunit,'(a)')'</DataArray>'
738  enddo
739  write(qunit,'(a)')'</PointData>'
740 
741  write(qunit,'(a)')'<Points>'
742  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
743  ! write cell corner coordinates in a backward dimensional loop, always 3D output
744  {do ix^db=ixcmin^db,ixcmax^db \}
745  x_vtk(1:3)=zero;
746  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
747  write(qunit,'(3(1pe14.6))') x_vtk
748  {end do \}
749  write(qunit,'(a)')'</DataArray>'
750  write(qunit,'(a)')'</Points>'
751  case('vtuCC')
752  ! we write out every grid as one VTK PIECE
753  write(qunit,'(a,i7,a,i7,a)') &
754  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
755  write(qunit,'(a)')'<CellData>'
756  do iw=1,nw+nwauxio
757  if(iw<=nw) then
758  if(.not.w_write(iw)) cycle
759  endif
760 
761  write(qunit,'(a,a,a)')&
762  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
763  write(qunit,'(200(1pe14.6))') {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
764  write(qunit,'(a)')'</DataArray>'
765  enddo
766  write(qunit,'(a)')'</CellData>'
767 
768  write(qunit,'(a)')'<Points>'
769  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
770  ! write cell corner coordinates in a backward dimensional loop, always 3D output
771  {do ix^db=ixcmin^db,ixcmax^db \}
772  x_vtk(1:3)=zero;
773  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
774  write(qunit,'(3(1pe14.6))') x_vtk
775  {end do \}
776  write(qunit,'(a)')'</DataArray>'
777  write(qunit,'(a)')'</Points>'
778  end select
779 
780 
781  write(qunit,'(a)')'<Cells>'
782 
783  ! connectivity part
784  write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
785  call save_connvtk(qunit,igrid)
786  write(qunit,'(a)')'</DataArray>'
787 
788  ! offsets data array
789  write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
790  do icel=1,nc
791  write(qunit,'(i7)') icel*(2**^nd)
792  end do
793  write(qunit,'(a)')'</DataArray>'
794 
795  ! VTK cell type data array
796  write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
797  ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
798  {^ifoned vtk_type=3 \}
799  {^iftwod vtk_type=8 \}
800  {^ifthreed vtk_type=11 \}
801  do icel=1,nc
802  write(qunit,'(i2)') vtk_type
803  enddo
804  write(qunit,'(a)')'</DataArray>'
805 
806  write(qunit,'(a)')'</Cells>'
807 
808  write(qunit,'(a)')'</Piece>'
809  endif
810  enddo
811  endif
812 enddo
813 
814 write(qunit,'(a)')'</UnstructuredGrid>'
815 write(qunit,'(a)')'</VTKFile>'
816 close(qunit)
817 
818 end subroutine unstructuredvtk
819 !====================================================================================
820 subroutine unstructuredvtkb(qunit)
822 ! output for vtu format to paraview, binary version output
823 ! not parallel, uses calc_grid to compute nwauxio variables
824 ! allows renormalizing using convert factors
828 
829 integer, intent(in) :: qunit
830 
831 double precision :: x_VTK(1:3)
832 
833 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
834 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
835 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
836 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
837 
838 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio):: wC_TMP
839 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
840 double precision :: normconv(0:nw+nwauxio)
841 
842 integer, allocatable :: intstatus(:,:)
843 integer :: ipe,igrid,level,icel,ixC^L,ixCC^L,Morton_no,Morton_length
844 integer :: nx^D,nxC^D,nc,np,VTK_type,ix^D,filenr
845 integer*8 :: offset
846 
847 integer:: k,iw
848 integer:: length,lengthcc,length_coords,length_conn,length_offsets
849 character:: buf
850 character(len=80):: filename
851 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
852 character(len=1024) :: outfilehead
853 
854 logical :: fileopen,cell_corner=.true.
855 logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
856 !-----------------------------------------------------------------------------
857 
858 normconv=one
859 morton_length=morton_stop(npe-1)-morton_start(0)+1
860 allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
861 allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
862 morton_aim=.false.
863 morton_aim_p=.false.
864 do morton_no=morton_start(mype),morton_stop(mype)
865  igrid=sfc_to_igrid(morton_no)
866  level=node(plevel_,igrid)
867  ! we can clip parts of the grid away, select variables, levels etc.
868  if(writelevel(level)) then
869  ! only output a grid when fully within clipped region selected
870  ! by writespshift array
871  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
872  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
873  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
874  morton_aim_p(morton_no)=.true.
875  end if
876  end if
877 end do
878 call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
879  icomm,ierrmpi)
880 select case(convert_type)
881  case('vtuB','vtuBmpi')
882  cell_corner=.true.
883  case('vtuBCC','vtuBCCmpi')
884  cell_corner=.false.
885 end select
886 if (mype /= 0) then
887  do morton_no=morton_start(mype),morton_stop(mype)
888  if(.not. morton_aim(morton_no)) cycle
889  igrid=sfc_to_igrid(morton_no)
890  call calc_x(igrid,xc,xcc)
891  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
892  ixc^l,ixcc^l,.true.)
893  itag=morton_no
894  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
895  if(cell_corner) then
896  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
897  else
898  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
899  endif
900  end do
901 
902 else
903  ! mype==0
904  offset=0
905 
906  inquire(qunit,opened=fileopen)
907  if(.not.fileopen)then
908  ! generate filename
909  filenr=snapshotini
910  if (autoconvert) filenr=snapshotnext
911  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
912  ! Open the file for the header part
913  open(qunit,file=filename,status='replace')
914  endif
915 
916  call getheadernames(wnamei,xandwnamei,outfilehead)
917 
918  ! generate xml header
919  write(qunit,'(a)')'<?xml version="1.0"?>'
920  write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
921  write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
922  write(qunit,'(a)')'<UnstructuredGrid>'
923  write(qunit,'(a)')'<FieldData>'
924  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
925  'NumberOfTuples="1" format="ascii">'
926  write(qunit,*) real(global_time*time_convert_factor)
927  write(qunit,'(a)')'</DataArray>'
928  write(qunit,'(a)')'</FieldData>'
929 
930  ! number of cells, number of corner points, per grid.
931  nx^d=ixmhi^d-ixmlo^d+1;
932  nxc^d=nx^d+1;
933  nc={nx^d*}
934  np={nxc^d*}
935 
936  length=np*size_real
937  lengthcc=nc*size_real
938 
939  length_coords=3*length
940  length_conn=2**^nd*size_int*nc
941  length_offsets=nc*size_int
942 
943  ! Note: using the w_write, writelevel, writespshift
944  do morton_no=morton_start(0),morton_stop(0)
945  if(.not. morton_aim(morton_no)) cycle
946  if(cell_corner) then
947  ! we write out every grid as one VTK PIECE
948  write(qunit,'(a,i7,a,i7,a)') &
949  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
950  write(qunit,'(a)')'<PointData>'
951  do iw=1,nw+nwauxio
952  if(iw<=nw) then
953  if(.not.w_write(iw)) cycle
954  endif
955 
956  write(qunit,'(a,a,a,i16,a)')&
957  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
958  '" format="appended" offset="',offset,'">'
959  write(qunit,'(a)')'</DataArray>'
960  offset=offset+length+size_int
961  enddo
962  write(qunit,'(a)')'</PointData>'
963 
964  write(qunit,'(a)')'<Points>'
965  write(qunit,'(a,i16,a)') &
966  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
967  ! write cell corner coordinates in a backward dimensional loop, always 3D output
968  offset=offset+length_coords+size_int
969  write(qunit,'(a)')'</Points>'
970  else
971  ! we write out every grid as one VTK PIECE
972  write(qunit,'(a,i7,a,i7,a)') &
973  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
974  write(qunit,'(a)')'<CellData>'
975  do iw=1,nw+nwauxio
976  if(iw<=nw) then
977  if(.not.w_write(iw)) cycle
978  endif
979 
980  write(qunit,'(a,a,a,i16,a)')&
981  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
982  '" format="appended" offset="',offset,'">'
983  write(qunit,'(a)')'</DataArray>'
984  offset=offset+lengthcc+size_int
985  enddo
986  write(qunit,'(a)')'</CellData>'
987 
988  write(qunit,'(a)')'<Points>'
989  write(qunit,'(a,i16,a)') &
990  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
991  ! write cell corner coordinates in a backward dimensional loop, always 3D output
992  offset=offset+length_coords+size_int
993  write(qunit,'(a)')'</Points>'
994  end if
995 
996  write(qunit,'(a)')'<Cells>'
997 
998  ! connectivity part
999  write(qunit,'(a,i16,a)')&
1000  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
1001  offset=offset+length_conn+size_int
1002 
1003  ! offsets data array
1004  write(qunit,'(a,i16,a)') &
1005  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
1006  offset=offset+length_offsets+size_int
1007 
1008  ! VTK cell type data array
1009  write(qunit,'(a,i16,a)') &
1010  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
1011  offset=offset+size_int+nc*size_int
1012 
1013  write(qunit,'(a)')'</Cells>'
1014 
1015  write(qunit,'(a)')'</Piece>'
1016  end do
1017  ! write metadata communicated from other processors
1018  if(npe>1)then
1019  do ipe=1, npe-1
1020  do morton_no=morton_start(ipe),morton_stop(ipe)
1021  if(.not. morton_aim(morton_no)) cycle
1022  if(cell_corner) then
1023  ! we write out every grid as one VTK PIECE
1024  write(qunit,'(a,i7,a,i7,a)') &
1025  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1026  write(qunit,'(a)')'<PointData>'
1027  do iw=1,nw+nwauxio
1028  if(iw<=nw) then
1029  if(.not.w_write(iw)) cycle
1030  endif
1031 
1032  write(qunit,'(a,a,a,i16,a)')&
1033  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
1034  '" format="appended" offset="',offset,'">'
1035  write(qunit,'(a)')'</DataArray>'
1036  offset=offset+length+size_int
1037  enddo
1038  write(qunit,'(a)')'</PointData>'
1039 
1040  write(qunit,'(a)')'<Points>'
1041  write(qunit,'(a,i16,a)') &
1042  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1043  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1044  offset=offset+length_coords+size_int
1045  write(qunit,'(a)')'</Points>'
1046  else
1047  ! we write out every grid as one VTK PIECE
1048  write(qunit,'(a,i7,a,i7,a)') &
1049  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1050  write(qunit,'(a)')'<CellData>'
1051  do iw=1,nw+nwauxio
1052  if(iw<=nw) then
1053  if(.not.w_write(iw)) cycle
1054  endif
1055 
1056  write(qunit,'(a,a,a,i16,a)')&
1057  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
1058  '" format="appended" offset="',offset,'">'
1059  write(qunit,'(a)')'</DataArray>'
1060  offset=offset+lengthcc+size_int
1061  enddo
1062  write(qunit,'(a)')'</CellData>'
1063 
1064  write(qunit,'(a)')'<Points>'
1065  write(qunit,'(a,i16,a)') &
1066  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1067  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1068  offset=offset+length_coords+size_int
1069  write(qunit,'(a)')'</Points>'
1070  end if
1071 
1072  write(qunit,'(a)')'<Cells>'
1073 
1074  ! connectivity part
1075  write(qunit,'(a,i16,a)')&
1076  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
1077  offset=offset+length_conn+size_int
1078 
1079  ! offsets data array
1080  write(qunit,'(a,i16,a)') &
1081  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
1082  offset=offset+length_offsets+size_int
1083 
1084  ! VTK cell type data array
1085  write(qunit,'(a,i16,a)') &
1086  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
1087  offset=offset+size_int+nc*size_int
1088 
1089  write(qunit,'(a)')'</Cells>'
1090 
1091  write(qunit,'(a)')'</Piece>'
1092  end do
1093  end do
1094  end if
1095 
1096  write(qunit,'(a)')'</UnstructuredGrid>'
1097  write(qunit,'(a)')'<AppendedData encoding="raw">'
1098  close(qunit)
1099  open(qunit,file=filename,access='stream',form='unformatted',position='append')
1100  buf='_'
1101  write(qunit) trim(buf)
1102 
1103  do morton_no=morton_start(0),morton_stop(0)
1104  if(.not. morton_aim(morton_no)) cycle
1105  igrid=sfc_to_igrid(morton_no)
1106  call calc_x(igrid,xc,xcc)
1107  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1108  ixc^l,ixcc^l,.true.)
1109  do iw=1,nw+nwauxio
1110  if(iw<=nw) then
1111  if(.not.w_write(iw)) cycle
1112  endif
1113  if(cell_corner) then
1114  write(qunit) length
1115  write(qunit) {(|}real(wC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixcmin^d,ixcmax^d)}
1116  else
1117  write(qunit) lengthcc
1118  write(qunit) {(|}real(wCC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixccmin^d,ixccmax^d)}
1119  end if
1120  enddo
1121 
1122  write(qunit) length_coords
1123  {do ix^db=ixcmin^db,ixcmax^db \}
1124  x_vtk(1:3)=zero;
1125  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1126  do k=1,3
1127  write(qunit) real(x_vtk(k))
1128  end do
1129  {end do \}
1130 
1131  write(qunit) length_conn
1132  {do ix^db=1,nx^db\}
1133  {^ifoned write(qunit)ix1-1,ix1 \}
1134  {^iftwod
1135  write(qunit)(ix2-1)*nxc1+ix1-1, &
1136  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1137  \}
1138  {^ifthreed
1139  write(qunit)&
1140  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1141  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1142  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1143  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1144  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1145  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1146  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1147  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1148  \}
1149  {end do\}
1150 
1151  write(qunit) length_offsets
1152  do icel=1,nc
1153  write(qunit) icel*(2**^nd)
1154  end do
1155 
1156 
1157  {^ifoned vtk_type=3 \}
1158  {^iftwod vtk_type=8 \}
1159  {^ifthreed vtk_type=11 \}
1160  write(qunit) size_int*nc
1161  do icel=1,nc
1162  write(qunit) vtk_type
1163  end do
1164  end do
1165  allocate(intstatus(mpi_status_size,1))
1166  if(npe>1)then
1167  ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d;
1168  ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
1169  do ipe=1, npe-1
1170  do morton_no=morton_start(ipe),morton_stop(ipe)
1171  if(.not. morton_aim(morton_no)) cycle
1172  itag=morton_no
1173  call mpi_recv(xc_tmp,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1174  if(cell_corner) then
1175  call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1176  else
1177  call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1178  end if
1179  do iw=1,nw+nwauxio
1180  if(iw<=nw) then
1181  if(.not.w_write(iw)) cycle
1182  endif
1183  if(cell_corner) then
1184  write(qunit) length
1185  write(qunit) {(|}real(wC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixcmin^d,ixcmax^d)}
1186  else
1187  write(qunit) lengthcc
1188  write(qunit) {(|}real(wCC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixccmin^d,ixccmax^d)}
1189  end if
1190  enddo
1191 
1192  write(qunit) length_coords
1193  {do ix^db=ixcmin^db,ixcmax^db \}
1194  x_vtk(1:3)=zero;
1195  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1196  do k=1,3
1197  write(qunit) real(x_vtk(k))
1198  end do
1199  {end do \}
1200 
1201  write(qunit) length_conn
1202  {do ix^db=1,nx^db\}
1203  {^ifoned write(qunit)ix1-1,ix1 \}
1204  {^iftwod
1205  write(qunit)(ix2-1)*nxc1+ix1-1, &
1206  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1207  \}
1208  {^ifthreed
1209  write(qunit)&
1210  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1211  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1212  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1213  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1214  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1215  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1216  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1217  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1218  \}
1219  {end do\}
1220 
1221  write(qunit) length_offsets
1222  do icel=1,nc
1223  write(qunit) icel*(2**^nd)
1224  end do
1225  {^ifoned vtk_type=3 \}
1226  {^iftwod vtk_type=8 \}
1227  {^ifthreed vtk_type=11 \}
1228  write(qunit) size_int*nc
1229  do icel=1,nc
1230  write(qunit) vtk_type
1231  end do
1232  end do
1233  end do
1234  end if
1235  close(qunit)
1236  open(qunit,file=filename,status='unknown',form='formatted',position='append')
1237  write(qunit,'(a)')'</AppendedData>'
1238  write(qunit,'(a)')'</VTKFile>'
1239  close(qunit)
1240  deallocate(intstatus)
1241 end if
1242 
1243 deallocate(morton_aim,morton_aim_p)
1244 if (npe>1) then
1245  call mpi_barrier(icomm,ierrmpi)
1246 endif
1247 
1248 end subroutine unstructuredvtkb
1249 !====================================================================================
1250 subroutine save_connvtk(qunit,igrid)
1252 ! this saves the basic line, pixel and voxel connectivity,
1253 ! as used by VTK file outputs for unstructured grid
1254 
1256 
1257 integer, intent(in) :: qunit, igrid
1258 
1259 integer :: nx^D, nxC^D, ix^D
1260 !-----------------------------------------------------------------------------
1261 nx^d=ixmhi^d-ixmlo^d+1;
1262 nxc^d=nx^d+1;
1263 
1264 {do ix^db=1,nx^db\}
1265  {^ifoned write(qunit,'(2(i7,1x))')ix1-1,ix1 \}
1266  {^iftwod
1267  write(qunit,'(4(i7,1x))')(ix2-1)*nxc1+ix1-1, &
1268  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1269  \}
1270  {^ifthreed
1271  write(qunit,'(8(i7,1x))')&
1272  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1273  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1274  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1275  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1276  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1277  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1278  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1279  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1280  \}
1281 {end do\}
1282 
1283 end subroutine save_connvtk
1284 !=============================================================================
1285 subroutine imagedatavtk_mpi(qunit)
1287 ! output for vti format to paraview, non-binary version output
1288 ! parallel, uses calc_grid to compute nwauxio variables
1289 ! allows renormalizing using convert factors
1290 ! allows skipping of w_write selected variables
1291 
1292 ! implementation such that length of ASCII output is identical when
1293 ! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
1294 
1297 use mod_calculate_xw
1298 
1299 integer, intent(in) :: qunit
1300 
1301 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
1302 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
1303 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1304 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1305 
1306 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
1307 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
1308 double precision, dimension(0:nw+nwauxio) :: normconv
1309 integer:: igrid,iigrid,level,ixC^L,ixCC^L
1310 integer:: NumGridsOnLevel(1:nlevelshi)
1311 integer :: nx^D
1312 
1313 character(len=80):: filename
1314 integer :: filenr
1315 
1316 integer, allocatable :: intstatus(:,:)
1317 logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
1318 
1319 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1320 character(len=1024) :: outfilehead
1321 
1322 logical :: fileopen
1323 integer :: ipe,Morton_no,Morton_length
1324 integer :: ixrvC^L, ixrvCC^L, siz_ind, ind_send(5*^nd), ind_recv(5*^nd)
1325 double precision :: origin(1:3), spacing(1:3)
1326 integer :: wholeExtent(1:6), ig^D
1327 type(tree_node_ptr) :: tree
1328 !-----------------------------------------------------------------------------
1329 if(levmin/=levmax) call mpistop('ImageData can only be used when levmin=levmax')
1330 
1331 normconv(0) = length_convert_factor
1332 normconv(1:nw) = w_convert_factor
1333 siz_ind=5*^nd
1334 morton_length=morton_stop(npe-1)-morton_start(0)+1
1335 allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
1336 allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
1337 morton_aim=.false.
1338 morton_aim_p=.false.
1339 do morton_no=morton_start(mype),morton_stop(mype)
1340  igrid=sfc_to_igrid(morton_no)
1341  level=node(plevel_,igrid)
1342  ! we can clip parts of the grid away, select variables, levels etc.
1343  if(writelevel(level)) then
1344  ! only output a grid when fully within clipped region selected
1345  ! by writespshift array
1346  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1347  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1348  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
1349  morton_aim_p(morton_no)=.true.
1350  end if
1351  end if
1352 end do
1353 call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
1354  icomm,ierrmpi)
1355 
1356 
1357 if (mype /= 0) then
1358  do morton_no=morton_start(mype),morton_stop(mype)
1359  if(.not. morton_aim(morton_no)) cycle
1360  igrid=sfc_to_igrid(morton_no)
1361  call calc_x(igrid,xc,xcc)
1362  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1363  ixc^l,ixcc^l,.true.)
1364  tree%node => igrid_to_node(igrid, mype)%node
1365  {^d& ig^d = tree%node%ig^d; }
1366  itag=morton_no
1367  ind_send=(/ ixc^l,ixcc^l, ig^d /)
1368  call mpi_send(ind_send,siz_ind,mpi_integer, 0,itag,icomm,ierrmpi)
1369  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
1370  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
1371  end do
1372 
1373 else
1374 
1375  inquire(qunit,opened=fileopen)
1376  if(.not.fileopen)then
1377  ! generate filename
1378  filenr=snapshotini
1379  if (autoconvert) filenr=snapshotnext
1380  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vti"
1381  ! Open the file for the header part
1382  open(qunit,file=filename,status='unknown',form='formatted')
1383  endif
1384 
1385 call getheadernames(wnamei,xandwnamei,outfilehead)
1386 
1387 ! number of cells per grid.
1388 nx^d=ixmhi^d-ixmlo^d+1;
1389 
1390 origin = 0
1391 {^d& origin(^d) = xprobmin^d*normconv(0); }
1392 spacing = zero
1393 {^d&spacing(^d) = dxlevel(^d)*normconv(0); }
1394 
1395 wholeextent = 0
1396 ! if we use writespshift, the whole extent has to be calculated:
1397 {^d&wholeextent(^d*2-1) = nx^d * ceiling(((xprobmax^d-xprobmin^d)*writespshift(^d,1)) &
1398  /(nx^d*dxlevel(^d))) \}
1399 {^d&wholeextent(^d*2) = nx^d * floor(((xprobmax^d-xprobmin^d)*(1.0d0-writespshift(^d,2))) &
1400  /(nx^d*dxlevel(^d))) \}
1401 
1402 ! generate xml header
1403 write(qunit,'(a)')'<?xml version="1.0"?>'
1404 write(qunit,'(a)',advance='no') '<VTKFile type="ImageData"'
1405 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1406 write(qunit,'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')' <ImageData Origin="',&
1407  origin,'" WholeExtent="',wholeextent,'" Spacing="',spacing,'">'
1408  write(qunit,'(a)')'<FieldData>'
1409  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1410  'NumberOfTuples="1" format="ascii">'
1411  write(qunit,*) real(global_time*time_convert_factor)
1412  write(qunit,'(a)')'</DataArray>'
1413  write(qunit,'(a)')'</FieldData>'
1414 
1415 ! write the data from proc 0
1416 do morton_no=morton_start(0),morton_stop(0)
1417  if(.not. morton_aim(morton_no)) cycle
1418  igrid=sfc_to_igrid(morton_no)
1419  tree%node => igrid_to_node(igrid, 0)%node
1420  {^d& ig^d = tree%node%ig^d; }
1421  call calc_x(igrid,xc,xcc)
1422  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1423  ixc^l,ixcc^l,.true.)
1424  call write_vti(qunit,ixg^ll,ixc^l,ixcc^l,ig^d,&
1425  nx^d,normconv,wnamei,wc_tmp,wcc_tmp)
1426 end do
1427 
1428 if(npe>1)then
1429  allocate(intstatus(mpi_status_size,1))
1430  do ipe=1, npe-1
1431  do morton_no=morton_start(ipe),morton_stop(ipe)
1432  if(.not. morton_aim(morton_no)) cycle
1433  itag=morton_no
1434  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1435  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
1436  ixrvccmin^d=ind_recv(2*^nd+^d);ixrvccmax^d=ind_recv(3*^nd+^d);
1437  ig^d=ind_recv(4*^nd+^d);
1438  call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1439  call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1440  call write_vti(qunit,ixg^ll,ixrvc^l,ixrvcc^l,ig^d,&
1441  nx^d,normconv,wnamei,wc_tmp,wcc_tmp)
1442  end do
1443  end do
1444 end if
1445 
1446 write(qunit,'(a)')'</ImageData>'
1447 write(qunit,'(a)')'</VTKFile>'
1448 close(qunit)
1449 if(npe>1) deallocate(intstatus)
1450 endif
1451 
1452 deallocate(morton_aim,morton_aim_p)
1453 if (npe>1) then
1454  call mpi_barrier(icomm,ierrmpi)
1455 endif
1456 
1457 end subroutine imagedatavtk_mpi
1458 !============================================================================
1459 subroutine punstructuredvtk_mpi(qunit)
1461 ! Write one pvtu and vtu files for each processor
1462 ! Otherwise like unstructuredvtk_mpi
1463 
1466 use mod_calculate_xw
1467 
1468 integer, intent(in) :: qunit
1469 !
1470 double precision, dimension(0:nw+nwauxio) :: normconv
1471 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
1472 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
1473 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1474 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1475 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
1476 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
1477 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1478 character(len=1024) :: outfilehead
1479 integer :: nx^D,nxC^D,nc,np, igrid,ixC^L,ixCC^L,level,Morton_no
1480 character(len=80) :: pfilename
1481 integer :: filenr
1482 logical :: fileopen,conv_grid
1483 !----------------------------------------------------------------------------
1484 
1485 ! Write pvtu-file:
1486 if (mype==0) then
1487  call write_pvtu(qunit)
1488 endif
1489 ! Now write the Source files:
1490 
1491 inquire(qunit,opened=fileopen)
1492 if(.not.fileopen)then
1493  ! generate filename
1494  filenr=snapshotini
1495  if (autoconvert) filenr=snapshotnext
1496  ! Open the file for the header part
1497  write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename),filenr,"p",mype,".vtu"
1498  open(qunit,file=pfilename,status='unknown',form='formatted')
1499 endif
1500 ! generate xml header
1501 write(qunit,'(a)')'<?xml version="1.0"?>'
1502 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
1503 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1504 write(qunit,'(a)')' <UnstructuredGrid>'
1505 write(qunit,'(a)')'<FieldData>'
1506 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1507  'NumberOfTuples="1" format="ascii">'
1508 write(qunit,*) real(global_time*time_convert_factor)
1509 write(qunit,'(a)')'</DataArray>'
1510 write(qunit,'(a)')'</FieldData>'
1511 
1512 call getheadernames(wnamei,xandwnamei,outfilehead)
1513 
1514 ! number of cells, number of corner points, per grid.
1515 nx^d=ixmhi^d-ixmlo^d+1;
1516 nxc^d=nx^d+1;
1517 nc={nx^d*}
1518 np={nxc^d*}
1519 
1520 ! Note: using the w_write, writelevel, writespshift
1521 ! we can clip parts of the grid away, select variables, levels etc.
1522 do level=levmin,levmax
1523  if (.not.writelevel(level)) cycle
1524  do morton_no=morton_start(mype),morton_stop(mype)
1525  igrid=sfc_to_igrid(morton_no)
1526  if (node(plevel_,igrid)/=level) cycle
1527 
1528  ! only output a grid when fully within clipped region selected
1529  ! by writespshift array
1530  conv_grid=({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1531  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1532  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})
1533  if (.not.conv_grid) cycle
1534 
1535  call calc_x(igrid,xc,xcc)
1536  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1537  ixc^l,ixcc^l,.true.)
1538 
1539  call write_vtk(qunit,ixg^ll,ixc^l,ixcc^l,igrid,nc,np,nx^d,nxc^d,&
1540  normconv,wnamei,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp)
1541  end do ! Morton_no loop
1542 end do ! level loop
1543 
1544  write(qunit,'(a)')' </UnstructuredGrid>'
1545  write(qunit,'(a)')'</VTKFile>'
1546  close(qunit)
1547 
1548 if (npe>1) then
1549  call mpi_barrier(icomm,ierrmpi)
1550 endif
1551 end subroutine punstructuredvtk_mpi
1552 !============================================================================
1553 subroutine unstructuredvtk_mpi(qunit)
1555 ! output for vtu format to paraview, non-binary version output
1556 ! parallel, uses calc_grid to compute nwauxio variables
1557 ! allows renormalizing using convert factors
1558 ! allows skipping of w_write selected variables
1559 
1560 ! implementation such that length of ASCII output is identical when
1561 ! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
1562 
1565 use mod_calculate_xw
1566 
1567 integer, intent(in) :: qunit
1568 
1569 double precision :: x_VTK(1:3)
1570 
1571 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
1572 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
1573 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1574 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1575 
1576 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
1577 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
1578 double precision, dimension(0:nw+nwauxio) :: normconv
1579 integer:: igrid,iigrid,level,ixC^L,ixCC^L
1580 integer:: NumGridsOnLevel(1:nlevelshi)
1581 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,ix^D
1582 
1583 character(len=80):: filename
1584 integer :: filenr
1585 
1586 integer, allocatable :: intstatus(:,:)
1587 
1588 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1589 character(len=1024) :: outfilehead
1590 
1591 logical :: fileopen,conv_grid,cond_grid_recv
1592 integer :: ipe,Morton_no,siz_ind
1593 integer :: ind_send(4*^nd),ind_recv(4*^nd)
1594 integer :: levmin_recv,levmax_recv,level_recv,igrid_recv,ixrvC^L,ixrvCC^L
1595 !-----------------------------------------------------------------------------
1596 if (mype==0) then
1597  inquire(qunit,opened=fileopen)
1598  if(.not.fileopen)then
1599  ! generate filename
1600  filenr=snapshotini
1601  if (autoconvert) filenr=snapshotnext
1602  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
1603  ! Open the file for the header part
1604  open(qunit,file=filename,status='unknown',form='formatted')
1605  endif
1606  ! generate xml header
1607  write(qunit,'(a)')'<?xml version="1.0"?>'
1608  write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
1609  write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1610  write(qunit,'(a)')'<UnstructuredGrid>'
1611  write(qunit,'(a)')'<FieldData>'
1612  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1613  'NumberOfTuples="1" format="ascii">'
1614  write(qunit,*) real(global_time*time_convert_factor)
1615  write(qunit,'(a)')'</DataArray>'
1616  write(qunit,'(a)')'</FieldData>'
1617 end if
1618 
1619 call getheadernames(wnamei,xandwnamei,outfilehead)
1620 ! number of cells, number of corner points, per grid.
1621 nx^d=ixmhi^d-ixmlo^d+1;
1622 nxc^d=nx^d+1;
1623 nc={nx^d*}
1624 np={nxc^d*}
1625 
1626 ! all slave processors send their minmal/maximal levels
1627 if (mype/=0) then
1628  if (morton_stop(mype)==0) call mpistop("nultag")
1629  itag=1000*morton_stop(mype)
1630  !print *,'ype,itag for levmin=',mype,itag,levmin
1631  call mpi_send(levmin,1,mpi_integer, 0,itag,icomm,ierrmpi)
1632  itag=2000*morton_stop(mype)
1633  !print *,'mype,itag for levmax=',mype,itag,levmax
1634  call mpi_send(levmax,1,mpi_integer, 0,itag,icomm,ierrmpi)
1635 end if
1636 
1637 
1638 ! Note: using the w_write, writelevel, writespshift
1639 ! we can clip parts of the grid away, select variables, levels etc.
1640 do level=levmin,levmax
1641  if (.not.writelevel(level)) cycle
1642  do morton_no=morton_start(mype),morton_stop(mype)
1643  igrid=sfc_to_igrid(morton_no)
1644  if (mype/=0)then
1645  itag=morton_no
1646  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
1647  itag=igrid
1648  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
1649  end if
1650  if (node(plevel_,igrid)/=level) cycle
1651 
1652  ! only output a grid when fully within clipped region selected
1653  ! by writespshift array
1654  conv_grid=({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1655  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1656  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})
1657  if (mype/=0)then
1658  call mpi_send(conv_grid,1,mpi_logical,0,itag,icomm,ierrmpi)
1659  end if
1660  if (.not.conv_grid) cycle
1661 
1662  call calc_x(igrid,xc,xcc)
1663  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1664  ixc^l,ixcc^l,.true.)
1665 
1666  if (mype/=0) then
1667  itag=morton_no
1668  ind_send=(/ ixc^l,ixcc^l /)
1669  siz_ind=4*^nd
1670  call mpi_send(ind_send,siz_ind,mpi_integer, 0,itag,icomm,ierrmpi)
1671  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
1672 
1673  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
1674  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
1675  itag=igrid
1676  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
1677  call mpi_send(xcc_tmp,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
1678  else
1679  call write_vtk(qunit,ixg^ll,ixc^l,ixcc^l,igrid,nc,np,nx^d,nxc^d,&
1680  normconv,wnamei,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp)
1681  end if
1682  end do ! Morton_no loop
1683 end do ! level loop
1684 
1685 
1686 if (mype==0) then
1687  allocate(intstatus(mpi_status_size,1))
1688  if(npe>1)then
1689  do ipe=1,npe-1
1690  itag=1000*morton_stop(ipe)
1691  call mpi_recv(levmin_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1692  !!print *,'mype RECEIVES,itag for levmin=',mype,itag,levmin_recv
1693  itag=2000*morton_stop(ipe)
1694  call mpi_recv(levmax_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1695  !!print *,'mype RECEIVES itag for levmax=',mype,itag,levmax_recv
1696  do level=levmin_recv,levmax_recv
1697  if (.not.writelevel(level)) cycle
1698  do morton_no=morton_start(ipe),morton_stop(ipe)
1699  itag=morton_no
1700  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1701  itag=igrid_recv
1702  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1703  if (level_recv/=level) cycle
1704 
1705  call mpi_recv(cond_grid_recv,1,mpi_logical, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1706  if(.not.cond_grid_recv)cycle
1707 
1708  itag=morton_no
1709  siz_ind=4*^nd
1710  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1711  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
1712  ixrvccmin^d=ind_recv(2*^nd+^d);ixrvccmax^d=ind_recv(3*^nd+^d);
1713  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
1714  ,icomm,intstatus(:,1),ierrmpi)
1715 
1716  call mpi_recv(wc_tmp_recv,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1717  call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1718 
1719  itag=igrid_recv
1720  call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1721  call mpi_recv(xcc_tmp_recv,1,type_block_xcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1722  call write_vtk(qunit,ixg^ll,ixrvc^l,ixrvcc^l,igrid_recv,&
1723  nc,np,nx^d,nxc^d,normconv,wnamei,&
1724  xc_tmp_recv,xcc_tmp_recv,wc_tmp_recv,wcc_tmp_recv)
1725  enddo ! Morton_no loop
1726  enddo ! level loop
1727  enddo ! processor loop
1728  endif ! multiple processors
1729  write(qunit,'(a)')'</UnstructuredGrid>'
1730  write(qunit,'(a)')'</VTKFile>'
1731  close(qunit)
1732 endif
1733 
1734 if (npe>1) then
1735  call mpi_barrier(icomm,ierrmpi)
1736  if(mype==0)deallocate(intstatus)
1737 endif
1738 
1739 end subroutine unstructuredvtk_mpi
1740 !============================================================================
1741 subroutine write_vtk(qunit,ixI^L,ixC^L,ixCC^L,igrid,nc,np,nx^D,nxC^D,&
1742  normconv,wnamei,xC,xCC,wC,wCC)
1745 
1746 integer, intent(in) :: qunit
1747 integer, intent(in) :: ixI^L,ixC^L,ixCC^L
1748 integer, intent(in) :: igrid,nc,np,nx^D,nxC^D
1749 double precision, intent(in) :: normconv(0:nw+nwauxio)
1750 character(len=name_len), intent(in):: wnamei(1:nw+nwauxio)
1751 
1752 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1753 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1754 
1755 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC
1756 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC
1757 
1758 double precision :: x_VTK(1:3)
1759 integer :: iw,ix^D,icel,VTK_type
1760 !----------------------------------------------------------------------------
1761 
1762 select case(convert_type)
1763  case('vtumpi','pvtumpi')
1764  ! we write out every grid as one VTK PIECE
1765  write(qunit,'(a,i7,a,i7,a)') &
1766  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1767  write(qunit,'(a)')'<PointData>'
1768  do iw=1,nw+nwauxio
1769  if(iw<=nw) then
1770  if(.not.w_write(iw)) cycle
1771  endif
1772 
1773  write(qunit,'(a,a,a)')&
1774  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
1775  write(qunit,'(200(1pe14.6))') {(|}wc(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
1776  write(qunit,'(a)')'</DataArray>'
1777  enddo
1778  write(qunit,'(a)')'</PointData>'
1779 
1780  write(qunit,'(a)')'<Points>'
1781  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
1782  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1783  {do ix^db=ixcmin^db,ixcmax^db \}
1784  x_vtk(1:3)=zero;
1785  x_vtk(1:ndim)=xc(ix^d,1:ndim)*normconv(0);
1786  write(qunit,'(3(1pe14.6))') x_vtk
1787  {end do \}
1788  write(qunit,'(a)')'</DataArray>'
1789  write(qunit,'(a)')'</Points>'
1790 
1791  case('vtuCCmpi','pvtuCCmpi')
1792  ! we write out every grid as one VTK PIECE
1793  write(qunit,'(a,i7,a,i7,a)') &
1794  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1795  write(qunit,'(a)')'<CellData>'
1796  do iw=1,nw+nwauxio
1797  if(iw<=nw) then
1798  if(.not.w_write(iw)) cycle
1799  endif
1800  write(qunit,'(a,a,a)')&
1801  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
1802  write(qunit,'(200(1pe14.6))') {(|}wcc(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
1803  write(qunit,'(a)')'</DataArray>'
1804  enddo
1805  write(qunit,'(a)')'</CellData>'
1806 
1807  write(qunit,'(a)')'<Points>'
1808  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
1809  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1810  {do ix^db=ixcmin^db,ixcmax^db \}
1811  x_vtk(1:3)=zero;
1812  x_vtk(1:ndim)=xc(ix^d,1:ndim)*normconv(0);
1813  write(qunit,'(3(1pe14.6))') x_vtk
1814  {end do \}
1815  write(qunit,'(a)')'</DataArray>'
1816  write(qunit,'(a)')'</Points>'
1817 end select
1818 
1819 write(qunit,'(a)')'<Cells>'
1820 
1821 ! connectivity part
1822 write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
1823 call save_connvtk(qunit,igrid)
1824 write(qunit,'(a)')'</DataArray>'
1825 
1826 ! offsets data array
1827 write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
1828 do icel=1,nc
1829  write(qunit,'(i7)') icel*(2**^nd)
1830 end do
1831 write(qunit,'(a)')'</DataArray>'
1832 
1833 ! VTK cell type data array
1834 write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
1835 ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
1836 {^ifoned vtk_type=3 \}
1837 {^iftwod vtk_type=8 \}
1838 {^ifthreed vtk_type=11 \}
1839 do icel=1,nc
1840  write(qunit,'(i2)') vtk_type
1841 enddo
1842 write(qunit,'(a)')'</DataArray>'
1843 
1844 write(qunit,'(a)')'</Cells>'
1845 
1846 write(qunit,'(a)')'</Piece>'
1847 
1848 end subroutine write_vtk
1849 !============================================================================
1850 subroutine write_vti(qunit,ixI^L,ixC^L,ixCC^L,ig^D,nx^D,&
1851  normconv,wnamei,wC,wCC)
1853 
1854 integer, intent(in) :: qunit
1855 integer, intent(in) :: ixI^L,ixC^L,ixCC^L
1856 integer, intent(in) :: ig^D,nx^D
1857 double precision, intent(in) :: normconv(0:nw+nwauxio)
1858 character(len=name_len), intent(in):: wnamei(1:nw+nwauxio)
1859 
1860 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC
1861 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC
1862 
1863 integer :: iw,ix^D
1864 integer :: extent(1:6)
1865 !----------------------------------------------------------------------------
1866 
1867 extent = 0
1868 {^d& extent(^d*2-1) = (ig^d-1) * nx^d; }
1869 {^d& extent(^d*2) = (ig^d) * nx^d; }
1870 
1871 
1872 select case(convert_type)
1873  case('vtimpi','pvtimpi')
1874  ! we write out every grid as one VTK PIECE
1875  write(qunit,'(a,6(i10),a)') &
1876  '<Piece Extent="',extent,'">'
1877  write(qunit,'(a)')'<PointData>'
1878  do iw=1,nw+nwauxio
1879  if(iw<=nw) then
1880  if(.not.w_write(iw)) cycle
1881  endif
1882 
1883  write(qunit,'(a,a,a)')&
1884  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
1885  write(qunit,'(200(1pe20.12))') {(|}wc(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
1886  write(qunit,'(a)')'</DataArray>'
1887  enddo
1888  write(qunit,'(a)')'</PointData>'
1889 
1890  case('vtiCCmpi','pvtiCCmpi')
1891  ! we write out every grid as one VTK PIECE
1892  write(qunit,'(a,6(i10),a)') &
1893  '<Piece Extent="',extent,'">'
1894  write(qunit,'(a)')'<CellData>'
1895  do iw=1,nw+nwauxio
1896  if(iw<=nw) then
1897  if(.not.w_write(iw)) cycle
1898  endif
1899  write(qunit,'(a,a,a)')&
1900  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
1901  write(qunit,'(200(1pe20.12))') {(|}wcc(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
1902  write(qunit,'(a)')'</DataArray>'
1903  enddo
1904  write(qunit,'(a)')'</CellData>'
1905 end select
1906 
1907 write(qunit,'(a)')'</Piece>'
1908 
1909 end subroutine write_vti
1910 !=============================================================================
1911 subroutine write_pvtu(qunit)
1914 use mod_calculate_xw
1915 
1916 integer, intent(in) :: qunit
1917 
1918 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio),outtype
1919 character(len=1024) :: outfilehead
1920 character(len=80) :: filename,pfilename
1921 integer :: filenr,iw,ipe,iscalars
1922 logical :: fileopen
1923 
1924 select case(convert_type)
1925 case('pvtumpi','pvtuBmpi')
1926  outtype="PPointData"
1927 case('pvtuCCmpi','pvtuBCCmpi')
1928  outtype="PCellData"
1929 end select
1930 inquire(qunit,opened=fileopen)
1931 if(.not.fileopen)then
1932  ! generate filename
1933  filenr=snapshotini
1934  if (autoconvert) filenr=snapshotnext
1935  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".pvtu"
1936  ! Open the file
1937  open(qunit,file=filename,status='unknown',form='formatted')
1938 endif
1939 
1940 call getheadernames(wnamei,xandwnamei,outfilehead)
1941 
1942 ! Get the default selection:
1943 iscalars=1
1944 do iw=nw,1, -1
1945  if (w_write(iw)) iscalars=iw
1946 end do
1947 
1948 
1949 ! generate xml header
1950 write(qunit,'(a)')'<?xml version="1.0"?>'
1951 write(qunit,'(a)',advance='no') '<VTKFile type="PUnstructuredGrid"'
1952 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1953 write(qunit,'(a)')' <PUnstructuredGrid GhostLevel="0">'
1954 ! Either celldata or pointdata:
1955 write(qunit,'(a,a,a,a,a)')&
1956  ' <',trim(outtype),' Scalars="',trim(wnamei(iscalars))//'">'
1957 do iw=1,nw
1958  if(.not.w_write(iw))cycle
1959  write(qunit,'(a,a,a)')&
1960  ' <PDataArray type="Float32" Name="',trim(wnamei(iw)),'"/>'
1961 end do
1962 do iw=nw+1,nw+nwauxio
1963  write(qunit,'(a,a,a)')&
1964  ' <PDataArray type="Float32" Name="',trim(wnamei(iw)),'"/>'
1965 end do
1966 write(qunit,'(a,a,a)')' </',trim(outtype),'>'
1967 
1968 write(qunit,'(a)')' <PPoints>'
1969 write(qunit,'(a)')' <PDataArray type="Float32" NumberOfComponents="3"/>'
1970 write(qunit,'(a)')' </PPoints>'
1971 
1972 do ipe=0,npe-1
1973  write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename(&
1974  index(base_filename, '/', back = .true.)+1:&
1975  len(base_filename))),filenr,"p",&
1976  ipe,".vtu"
1977  write(qunit,'(a,a,a)')' <Piece Source="',trim(pfilename),'"/>'
1978 end do
1979 write(qunit,'(a)')' </PUnstructuredGrid>'
1980 write(qunit,'(a)')'</VTKFile>'
1981 
1982 close(qunit)
1983 
1984 end subroutine write_pvtu
1985 !=============================================================================
1986 subroutine tecplot_mpi(qunit)
1988 ! output for tecplot (ASCII format)
1989 ! parallel, uses calc_grid to compute nwauxio variables
1990 ! allows renormalizing using convert factors
1991 
1992 ! the current implementation is such that tecplotmpi and tecplotCCmpi will
1993 ! create different length output ASCII files when used on 1 versus multiple CPUs
1994 ! in fact, on 1 CPU, there will be as many zones as there are levels
1995 ! on multiple CPUs, there will be a number of zones up to the number of
1996 ! levels times the number of CPUs (can be less, when some level not on a CPU)
1997 
2000 use mod_calculate_xw
2001 
2002 integer, intent(in) :: qunit
2003 
2004 integer:: igrid,iigrid,level,igonlevel,iw,idim,ix^D
2005 integer:: NumGridsOnLevel(1:nlevelshi)
2006 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,ixC^L,ixCC^L
2007 integer :: nodesonlevelmype,elemsonlevelmype
2008 
2009 integer :: nodes, elems
2010 
2011 integer, allocatable :: intstatus(:,:)
2012 
2013 double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
2014 
2015 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
2016 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
2017 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
2018 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
2019 
2020 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
2021 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
2022 double precision, dimension(0:nw+nwauxio) :: normconv
2023 logical :: fileopen,first
2024 integer :: Morton_no,ipe,levmin_recv,levmax_recv,igrid_recv,level_recv
2025 integer :: ixrvC^L,ixrvCC^L
2026 integer :: ind_send(2*^nd),ind_recv(2*^nd),siz_ind,igonlevel_recv
2027 integer :: NumGridsOnLevel_mype(1:nlevelshi,0:npe-1)
2028 character(len=80) :: filename
2029 integer :: filenr
2030 character(len=1024) :: tecplothead
2031 
2032 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
2033 character(len=1024) :: outfilehead
2034 !-----------------------------------------------------------------------------
2035 if(nw/=count(w_write(1:nw)))then
2036  if(mype==0) print *,'tecplot_mpi does not use w_write=F'
2037  call mpistop('w_write, tecplot')
2038 end if
2039 
2040 if(nocartesian)then
2041  if(mype==0) print *,'tecplot_mpi with nocartesian'
2042 endif
2043 
2044 master_cpu_open : if (mype == 0) then
2045  inquire(qunit,opened=fileopen)
2046  if (.not.fileopen) then
2047  ! generate filename
2048  filenr=snapshotini
2049  if (autoconvert) filenr=snapshotnext
2050  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".plt"
2051  open(qunit,file=filename,status='unknown')
2052  end if
2053 
2054  call getheadernames(wnamei,xandwnamei,outfilehead)
2055 
2056  write(tecplothead,'(a)') "VARIABLES = "//trim(outfilehead)
2057  write(qunit,'(a)') tecplothead(1:len_trim(tecplothead))
2058 end if master_cpu_open
2059 
2060 
2061 ! determine overall number of grids per level, and the same info per CPU
2062 numgridsonlevel(1:nlevelshi)=0
2063 do level=levmin,levmax
2064  numgridsonlevel(level)=0
2065  do morton_no=morton_start(mype),morton_stop(mype)
2066  igrid = sfc_to_igrid(morton_no)
2067  if (node(plevel_,igrid)/=level) cycle
2068  numgridsonlevel(level)=numgridsonlevel(level)+1
2069  end do
2070  numgridsonlevel_mype(level,0:npe-1)=0
2071  numgridsonlevel_mype(level,mype) = numgridsonlevel(level)
2072  call mpi_allreduce(mpi_in_place,numgridsonlevel_mype(level,0:npe-1),npe,mpi_integer,&
2073  mpi_max,icomm,ierrmpi)
2074  call mpi_allreduce(mpi_in_place,numgridsonlevel(level),1,mpi_integer,mpi_sum, &
2075  icomm,ierrmpi)
2076 end do
2077 
2078 
2079 !!do level=levmin,levmax
2080 !! print *,'mype, level en NumGridsOnLevel_mype(level,0:npe-1)=', &
2081 !! mype,level,NumGridsOnLevel_mype(level,0:npe-1)
2082 !! print *,'mype, level en NumGridsOnLevel(level)=', &
2083 !! mype,level,NumGridsOnLevel(level)
2084 !!enddo
2085 
2086 
2087 nx^d=ixmhi^d-ixmlo^d+1;
2088 nxc^d=nx^d+1;
2089 
2090 if(mype==0.and.npe>1) allocate(intstatus(mpi_status_size,1))
2091 
2092 {^ifoned
2093 if(convert_type=='teclinempi') then
2094  nodes=0
2095  elems=0
2096  do level=levmin,levmax
2097  nodes=nodes + numgridsonlevel(level)*{nxc^d*}
2098  elems=elems + numgridsonlevel(level)*{nx^d*}
2099  enddo
2100 
2101  if (mype==0) write(qunit,"(a,i7,a,1pe12.5,a)") &
2102  'ZONE T="all levels", I=',elems, &
2103  ', SOLUTIONTIME=',global_time*time_convert_factor,', F=POINT'
2104 
2105  igonlevel=0
2106  do morton_no=morton_start(mype),morton_stop(mype)
2107  igrid = sfc_to_igrid(morton_no)
2108  call calc_x(igrid,xc,xcc)
2109  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,ixc^l,ixcc^l,.true.)
2110  if (mype==0) then
2111  {do ix^db=ixccmin^db,ixccmax^db\}
2112  x_tec(1:ndim)=xcc_tmp(ix^d,1:ndim)*normconv(0)
2113  w_tec(1:nw+nwauxio)=wcc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2114  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2115  {end do\}
2116  else if (mype/=0) then
2117  itag=morton_no
2118  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2119  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision,0,itag,icomm,ierrmpi)
2120  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
2121  call mpi_send(xcc_tmp,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
2122  end if
2123  enddo
2124  if (mype==0) then
2125  do ipe=1,npe-1
2126  do morton_no=morton_start(ipe),morton_stop(ipe)
2127  itag=morton_no
2128  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2129  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,&
2130  itag,icomm,intstatus(:,1),ierrmpi)
2131  call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,&
2132  icomm,intstatus(:,1),ierrmpi)
2133  call mpi_recv(xcc_tmp_recv,1,type_block_xcc_io, ipe,itag,&
2134  icomm,intstatus(:,1),ierrmpi)
2135  {do ix^db=ixccmin^db,ixccmax^db\}
2136  x_tec(1:ndim)=xcc_tmp_recv(ix^d,1:ndim)*normconv(0)
2137  w_tec(1:nw+nwauxio)=wcc_tmp_recv(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2138  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2139  {end do\}
2140  end do
2141  end do
2142  close(qunit)
2143  end if
2144 else
2145 }
2146 
2147 if (mype/=0) then
2148  itag=1000*morton_stop(mype)
2149  call mpi_send(levmin,1,mpi_integer, 0,itag,icomm,ierrmpi)
2150  itag=2000*morton_stop(mype)
2151  call mpi_send(levmax,1,mpi_integer, 0,itag,icomm,ierrmpi)
2152 end if
2153 
2154 do level=levmin,levmax
2155  nodesonlevelmype=numgridsonlevel_mype(level,mype)*{nxc^d*}
2156  elemsonlevelmype=numgridsonlevel_mype(level,mype)*{nx^d*}
2157  nodesonlevel=numgridsonlevel(level)*{nxc^d*}
2158  elemsonlevel=numgridsonlevel(level)*{nx^d*}
2159  ! for all tecplot variants coded up here, we let the TECPLOT ZONES coincide
2160  ! with the AMR grid LEVEL. Other options would be
2161  ! let each grid define a zone: inefficient for TECPLOT internal workings
2162  ! hence not implemented
2163  ! let entire octree define 1 zone: no difference in interpolation
2164  ! properties across TECPLOT zones detected as yet, hence not done
2165  select case(convert_type)
2166  case('tecplotmpi')
2167  ! in this option, we store the corner coordinates, as well as the corner
2168  ! values of all variables (obtained by averaging). This allows POINT packaging,
2169  ! and thus we can save full grid info by using one call to calc_grid
2170  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2171  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
2172  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2173  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
2174  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2175  do morton_no=morton_start(mype),morton_stop(mype)
2176  igrid = sfc_to_igrid(morton_no)
2177  if (mype/=0)then
2178  itag=morton_no
2179  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2180  itag=igrid
2181  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2182  end if
2183  if (node(plevel_,igrid)/=level) cycle
2184  call calc_x(igrid,xc,xcc)
2185  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2186  ixc^l,ixcc^l,.true.)
2187  if (mype/=0) then
2188  itag=morton_no
2189  ind_send=(/ ixc^l /)
2190  siz_ind=2*^nd
2191  call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2192  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2193 
2194  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
2195  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
2196  else
2197  {do ix^db=ixcmin^db,ixcmax^db\}
2198  x_tec(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0)
2199  w_tec(1:nw+nwauxio)=wc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2200  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2201  {end do\}
2202  end if
2203  enddo
2204 
2205  case('tecplotCCmpi')
2206  ! in this option, we store the corner coordinates, and the cell center
2207  ! values of all variables. Due to this mix of corner/cell center, we must
2208  ! use BLOCK packaging, and thus we have enormous overhead by using
2209  ! calc_grid repeatedly to merely fill values of cell corner coordinates
2210  ! and cell center values per dimension, per variable
2211  if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
2212  if(nw+nwauxio==1)then
2213  ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
2214  ! and just set [ndim+1]
2215  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2216  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
2217  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2218  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2219  ndim+1,']=CELLCENTERED), ZONETYPE=', &
2220  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2221  else
2222  if(ndim+nw+nwauxio<10) then
2223  ! difference only in length of integer format specification for ndim+nw+nwauxio
2224  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2225  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
2226  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2227  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2228  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2229  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2230  else
2231  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2232  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
2233  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2234  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2235  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2236  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2237  endif
2238  endif
2239 
2240  do idim=1,ndim
2241  first=(idim==1)
2242  do morton_no=morton_start(mype),morton_stop(mype)
2243  igrid = sfc_to_igrid(morton_no)
2244  if (mype/=0)then
2245  itag=morton_no*idim
2246  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2247  itag=igrid*idim
2248  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2249  end if
2250  if (node(plevel_,igrid)/=level) cycle
2251 
2252  call calc_x(igrid,xc,xcc)
2253  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2254  ixc^l,ixcc^l,first)
2255  if (mype/=0)then
2256  ind_send=(/ ixc^l /)
2257  siz_ind=2*^nd
2258  itag=igrid*idim
2259  call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2260  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2261  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
2262  else
2263  write(qunit,fmt="(100(e14.6))") xc_tmp(ixc^s,idim)*normconv(0)
2264  end if
2265  enddo
2266  enddo
2267 
2268  do iw=1,nw+nwauxio
2269  do morton_no=morton_start(mype),morton_stop(mype)
2270  igrid = sfc_to_igrid(morton_no)
2271  if (mype/=0)then
2272  itag=morton_no*(ndim+iw)
2273  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2274  itag=igrid*(ndim+iw)
2275  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2276  end if
2277  if (node(plevel_,igrid)/=level) cycle
2278 
2279  call calc_x(igrid,xc,xcc)
2280  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2281  ixc^l,ixcc^l,.true.)
2282 
2283  if (mype/=0)then
2284  ind_send=(/ ixcc^l /)
2285  siz_ind=2*^nd
2286  itag=igrid*(ndim+iw)
2287  call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2288  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2289  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
2290  else
2291  write(qunit,fmt="(100(e14.6))") wcc_tmp(ixcc^s,iw)*normconv(iw)
2292  endif
2293  enddo
2294  enddo
2295  case default
2296  call mpistop('no such tecplot type')
2297  end select
2298 
2299 
2300  igonlevel=0
2301  do morton_no=morton_start(mype),morton_stop(mype)
2302  igrid = sfc_to_igrid(morton_no)
2303  if (mype/=0)then
2304  itag=morton_no
2305  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2306  itag=igrid
2307  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2308  end if
2309  if (node(plevel_,igrid)/=level) cycle
2310 
2311  igonlevel=igonlevel+1
2312  if (mype/=0)then
2313  itag=igrid
2314  call mpi_send(igonlevel,1,mpi_integer, 0,itag,icomm,ierrmpi)
2315  end if
2316  if(mype==0)then
2317  call save_conntec(qunit,igrid,igonlevel)
2318  endif
2319  end do
2320 end do
2321 
2322 if (mype==0) then
2323  if (npe>1) then
2324  do ipe=1,npe-1
2325  itag=1000*morton_stop(ipe)
2326  call mpi_recv(levmin_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2327  itag=2000*morton_stop(ipe)
2328  call mpi_recv(levmax_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2329  do level=levmin_recv,levmax_recv
2330  nodesonlevelmype=numgridsonlevel_mype(level,ipe)*{nxc^d*}
2331  elemsonlevelmype=numgridsonlevel_mype(level,ipe)*{nx^d*}
2332  nodesonlevel=numgridsonlevel(level)*{nxc^d*}
2333  elemsonlevel=numgridsonlevel(level)*{nx^d*}
2334  select case(convert_type)
2335  case('tecplotmpi')
2336  ! in this option, we store the corner coordinates, as well as the corner
2337  ! values of all variables (obtained by averaging). This allows POINT packaging,
2338  ! and thus we can save full grid info by using one call to calc_grid
2339  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2340  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
2341  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2342  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
2343  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2344  do morton_no=morton_start(ipe),morton_stop(ipe)
2345  itag=morton_no
2346  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2347  itag=igrid_recv
2348  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2349  if (level_recv/=level) cycle
2350 
2351  itag=morton_no
2352  siz_ind=2*^nd
2353  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,&
2354  icomm,intstatus(:,1),ierrmpi)
2355  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
2356  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2357  ,icomm,intstatus(:,1),ierrmpi)
2358 
2359  call mpi_recv(wc_tmp_recv,1,type_block_wc_io, ipe,itag,&
2360  icomm,intstatus(:,1),ierrmpi)
2361  call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,&
2362  icomm,intstatus(:,1),ierrmpi)
2363 
2364  {do ix^db=ixrvcmin^db,ixrvcmax^db\}
2365  x_tec(1:ndim)=xc_tmp_recv(ix^d,1:ndim)*normconv(0)
2366  w_tec(1:nw+nwauxio)=wc_tmp_recv(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2367  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2368  {end do\}
2369  end do
2370  case('tecplotCCmpi')
2371  ! in this option, we store the corner coordinates, and the cell center
2372  ! values of all variables. Due to this mix of corner/cell center, we must
2373  ! use BLOCK packaging, and thus we have enormous overhead by using
2374  ! calc_grid repeatedly to merely fill values of cell corner coordinates
2375  ! and cell center values per dimension, per variable
2376  if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
2377  if(nw+nwauxio==1)then
2378  ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
2379  ! and just set [ndim+1]
2380  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2381  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
2382  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2383  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2384  ndim+1,']=CELLCENTERED), ZONETYPE=', &
2385  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2386  else
2387  if(ndim+nw+nwauxio<10) then
2388  ! difference only in length of integer format specification for ndim+nw+nwauxio
2389  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2390  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
2391  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2392  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2393  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2394  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2395  else
2396  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2397  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
2398  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2399  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2400  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2401  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2402  endif
2403  endif
2404 
2405  do idim=1,ndim
2406  do morton_no=morton_start(ipe),morton_stop(ipe)
2407  itag=morton_no*idim
2408  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2409  itag=igrid_recv*idim
2410  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2411  if (level_recv/=level) cycle
2412 
2413  siz_ind=2*^nd
2414  itag=igrid_recv*idim
2415  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2416  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
2417  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2418  ,icomm,intstatus(:,1),ierrmpi)
2419  call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2420  write(qunit,fmt="(100(e14.6))") xc_tmp_recv(ixrvc^s,idim)*normconv(0)
2421  end do
2422  end do
2423 
2424  do iw=1,nw+nwauxio
2425  do morton_no=morton_start(ipe),morton_stop(ipe)
2426  itag=morton_no*(ndim+iw)
2427  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2428  itag=igrid_recv*(ndim+iw)
2429  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2430  if (level_recv/=level) cycle
2431 
2432  siz_ind=2*^nd
2433  itag=igrid_recv*(ndim+iw)
2434  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2435  ixrvccmin^d=ind_recv(^d);ixrvccmax^d=ind_recv(^nd+^d);
2436  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2437  ,icomm,intstatus(:,1),ierrmpi)
2438  call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2439  write(qunit,fmt="(100(e14.6))") wcc_tmp_recv(ixrvcc^s,iw)*normconv(iw)
2440  enddo
2441  enddo
2442  case default
2443  call mpistop('no such tecplot type')
2444  end select
2445 
2446  do morton_no=morton_start(ipe),morton_stop(ipe)
2447  itag=morton_no
2448  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2449  itag=igrid_recv
2450  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2451  if (level_recv/=level) cycle
2452 
2453  itag=igrid_recv
2454  call mpi_recv(igonlevel_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2455  call save_conntec(qunit,igrid_recv,igonlevel_recv)
2456  end do ! morton loop
2457  end do ! level loop
2458  end do ! ipe loop
2459  end if ! npe>1 if
2460 end if ! mype=0 if
2461 {^ifoned endif}
2462 
2463 if (npe>1) then
2464  call mpi_barrier(icomm,ierrmpi)
2465  if(mype==0)deallocate(intstatus)
2466 endif
2467 
2468 end subroutine tecplot_mpi
2469 !=============================================================================
2470 subroutine punstructuredvtkb_mpi(qunit)
2472 ! Write one pvtu and vtu files for each processor
2473 ! Otherwise like unstructuredvtk_mpi
2474 ! output for vtu format to paraview, binary version output
2475 ! uses calc_grid to compute nwauxio variables
2476 ! allows renormalizing using convert factors
2477 
2480 use mod_calculate_xw
2481 
2482 integer, intent(in) :: qunit
2483 
2484 double precision :: x_VTK(1:3)
2485 
2486 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
2487 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
2488 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
2489 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
2490 
2491 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
2492 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
2493 
2494 integer :: igrid,iigrid,level,igonlevel,icel,ixC^L,ixCC^L,Morton_no
2495 integer :: NumGridsOnLevel(1:nlevelshi)
2496 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,VTK_type,ix^D
2497 double precision :: normconv(0:nw+nwauxio)
2498 character(len=80) :: pfilename
2499 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
2500 character(len=1024) :: outfilehead
2501 
2502 integer*8 :: offset
2503 integer:: recsep,k,iw,filenr
2504 integer:: length,lengthcc,offset_points,offset_cells, &
2505  length_coords,length_conn,length_offsets
2506 character:: buf
2507 character(len=6):: bufform
2508 
2509 logical :: fileopen
2510 !-----------------------------------------------------------------------------
2511 
2512 ! Write pvtu-file:
2513 if (mype==0) then
2514  call write_pvtu(qunit)
2515 endif
2516 ! Now write the Source files:
2517 inquire(qunit,opened=fileopen)
2518 if(.not.fileopen)then
2519  ! generate filename
2520  filenr=snapshotnext-1
2521  if (autoconvert) filenr=snapshotnext
2522  ! Open the file for the header part
2523  write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename),filenr,"p",mype,".vtu"
2524  open(qunit,file=pfilename,status='unknown',form='formatted')
2525 endif
2526 ! generate xml header
2527 write(qunit,'(a)')'<?xml version="1.0"?>'
2528 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
2529 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
2530 write(qunit,'(a)')' <UnstructuredGrid>'
2531 write(qunit,'(a)')'<FieldData>'
2532 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
2533  'NumberOfTuples="1" format="ascii">'
2534 write(qunit,*) real(global_time*time_convert_factor)
2535 write(qunit,'(a)')'</DataArray>'
2536 write(qunit,'(a)')'</FieldData>'
2537 
2538 
2539 offset=0
2540 recsep=4
2541 
2542 call getheadernames(wnamei,xandwnamei,outfilehead)
2543 
2544 ! number of cells, number of corner points, per grid.
2545 nx^d=ixmhi^d-ixmlo^d+1;
2546 nxc^d=nx^d+1;
2547 nc={nx^d*}
2548 np={nxc^d*}
2549 
2550 length=np*size_real
2551 lengthcc=nc*size_real
2552 
2553 length_coords=3*length
2554 length_conn=2**^nd*size_int*nc
2555 length_offsets=nc*size_int
2556 
2557 ! Note: using the w_write, writelevel, writespshift
2558 ! we can clip parts of the grid away, select variables, levels etc.
2559 do level=levmin,levmax
2560  if (writelevel(level)) then
2561  do morton_no=morton_start(mype),morton_stop(mype)
2562  igrid=sfc_to_igrid(morton_no)
2563  if (node(plevel_,igrid)/=level) cycle
2564  ! only output a grid when fully within clipped region selected
2565  ! by writespshift array
2566  if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
2567  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
2568  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
2569  select case(convert_type)
2570  case('pvtuBmpi')
2571  ! we write out every grid as one VTK PIECE
2572  write(qunit,'(a,i7,a,i7,a)') &
2573  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2574  write(qunit,'(a)')'<PointData>'
2575  do iw=1,nw
2576  if(.not.w_write(iw))cycle
2577 
2578  write(qunit,'(a,a,a,i16,a)')&
2579  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2580  '" format="appended" offset="',offset,'">'
2581  write(qunit,'(a)')'</DataArray>'
2582  offset=offset+length+size_int
2583  enddo
2584  do iw=nw+1,nw+nwauxio
2585 
2586  write(qunit,'(a,a,a,i16,a)')&
2587  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2588  '" format="appended" offset="',offset,'">'
2589  write(qunit,'(a)')'</DataArray>'
2590  offset=offset+length+size_int
2591  enddo
2592  write(qunit,'(a)')'</PointData>'
2593 
2594  write(qunit,'(a)')'<Points>'
2595  write(qunit,'(a,i16,a)') &
2596  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
2597  ! write cell corner coordinates in a backward dimensional loop, always 3D output
2598  offset=offset+length_coords+size_int
2599  write(qunit,'(a)')'</Points>'
2600  case('pvtuBCCmpi')
2601  ! we write out every grid as one VTK PIECE
2602  write(qunit,'(a,i7,a,i7,a)') &
2603  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2604  write(qunit,'(a)')'<CellData>'
2605  do iw=1,nw
2606  if(.not.w_write(iw))cycle
2607 
2608  write(qunit,'(a,a,a,i16,a)')&
2609  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2610  '" format="appended" offset="',offset,'">'
2611  write(qunit,'(a)')'</DataArray>'
2612  offset=offset+lengthcc+size_int
2613  enddo
2614  do iw=nw+1,nw+nwauxio
2615 
2616  write(qunit,'(a,a,a,i16,a)')&
2617  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2618  '" format="appended" offset="',offset,'">'
2619  write(qunit,'(a)')'</DataArray>'
2620  offset=offset+lengthcc+size_int
2621  enddo
2622  write(qunit,'(a)')'</CellData>'
2623 
2624  write(qunit,'(a)')'<Points>'
2625  write(qunit,'(a,i16,a)') &
2626  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
2627  ! write cell corner coordinates in a backward dimensional loop, always 3D output
2628  offset=offset+length_coords+size_int
2629  write(qunit,'(a)')'</Points>'
2630  end select
2631 
2632 
2633  write(qunit,'(a)')'<Cells>'
2634 
2635  ! connectivity part
2636  write(qunit,'(a,i16,a)')&
2637  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
2638  offset=offset+length_conn+size_int
2639 
2640  ! offsets data array
2641  write(qunit,'(a,i16,a)') &
2642  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
2643  offset=offset+length_offsets+size_int
2644 
2645  ! VTK cell type data array
2646  write(qunit,'(a,i16,a)') &
2647  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
2648  offset=offset+size_int+nc*size_int
2649 
2650  write(qunit,'(a)')'</Cells>'
2651 
2652  write(qunit,'(a)')'</Piece>'
2653  endif
2654  enddo
2655  endif
2656 enddo
2657 
2658 write(qunit,'(a)')'</UnstructuredGrid>'
2659 write(qunit,'(a)')'<AppendedData encoding="raw">'
2660 
2661 close(qunit)
2662 ! next to make gfortran compiler happy, as it does not know
2663 ! form='binary' and produces error on compilation
2664 !bufform='binary'
2665 !open(qunit,file=pfilename,form=bufform,position='append')
2666 !This should in principle do also for gfortran (tested with gfortran 4.6.0 and Intel 11.1):
2667 open(qunit,file=pfilename,access='stream',form='unformatted',position='append')
2668 buf='_'
2669 write(qunit) trim(buf)
2670 
2671 do level=levmin,levmax
2672  if (writelevel(level)) then
2673  do morton_no=morton_start(mype),morton_stop(mype)
2674  igrid=sfc_to_igrid(morton_no)
2675  if (node(plevel_,igrid)/=level) cycle
2676  ! only output a grid when fully within clipped region selected
2677  ! by writespshift array
2678  if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
2679  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
2680  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
2681  call calc_x(igrid,xc,xcc)
2682  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2683  ixc^l,ixcc^l,.true.)
2684  do iw=1,nw
2685  if(.not.w_write(iw))cycle
2686  select case(convert_type)
2687  case('pvtuBmpi')
2688  write(qunit) length
2689  write(qunit) {(|}real(wC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixcmin^d,ixcmax^d)}
2690  case('pvtuBCCmpi')
2691  write(qunit) lengthcc
2692  write(qunit) {(|}real(wCC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixccmin^d,ixccmax^d)}
2693  end select
2694  enddo
2695  do iw=nw+1,nw+nwauxio
2696  select case(convert_type)
2697  case('pvtuBmpi')
2698  write(qunit) length
2699  write(qunit) {(|}real(wC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixcmin^d,ixcmax^d)}
2700  case('pvtuBCCmpi')
2701  write(qunit) lengthcc
2702  write(qunit) {(|}real(wCC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixccmin^d,ixccmax^d)}
2703  end select
2704  enddo
2705 
2706  write(qunit) length_coords
2707  {do ix^db=ixcmin^db,ixcmax^db \}
2708  x_vtk(1:3)=zero;
2709  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
2710  do k=1,3
2711  write(qunit) real(x_vtk(k))
2712  end do
2713  {end do \}
2714 
2715  write(qunit) length_conn
2716  {do ix^db=1,nx^db\}
2717  {^ifoned write(qunit)ix1-1,ix1 \}
2718  {^iftwod
2719  write(qunit)(ix2-1)*nxc1+ix1-1, &
2720  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
2721  \}
2722  {^ifthreed
2723  write(qunit)&
2724  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
2725  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
2726  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
2727  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
2728  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
2729  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
2730  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
2731  ix3*nxc2*nxc1+ ix2*nxc1+ix1
2732  \}
2733  {end do\}
2734 
2735  write(qunit) length_offsets
2736  do icel=1,nc
2737  write(qunit) icel*(2**^nd)
2738  end do
2739 
2740 
2741  {^ifoned vtk_type=3 \}
2742  {^iftwod vtk_type=8 \}
2743  {^ifthreed vtk_type=11 \}
2744  write(qunit) size_int*nc
2745  do icel=1,nc
2746  write(qunit) vtk_type
2747  enddo
2748  endif
2749  end do
2750  endif
2751 end do
2752 
2753 close(qunit)
2754 open(qunit,file=pfilename,status='unknown',form='formatted',position='append')
2755 
2756 write(qunit,'(a)')'</AppendedData>'
2757 write(qunit,'(a)')'</VTKFile>'
2758 close(qunit)
2759 
2760 end subroutine punstructuredvtkb_mpi
2761 !=============================================================================
subroutine write_pvtu(qunit)
Definition: convert.t:1912
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine get_euv_image(qunit)
subroutine imagedatavtk_mpi(qunit)
Definition: convert.t:1286
update ghost cells of all blocks including physical boundaries
logical, dimension(:), allocatable w_write
subroutine onegrid(qunit)
Definition: convert.t:254
logical saveprim
If true, convert from conservative to primitive variables in output.
integer max_blocks
The maximum number of grid blocks in a processor.
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 oneblock(qunit)
Definition: convert.t:58
integer type_block_wc_io
MPI type for IO: cell corner (wc) or cell center (wcc) variables.
subroutine getheadernames(wnamei, xandwnamei, outfilehead)
get all variables names
double precision, dimension(ndim, 2) writespshift
double precision time_convert_factor
Conversion factor for time unit.
character(len=std_len) convert_type
Which format to use when converting.
integer, parameter plevel_
subroutine write_vti(qunit, ixIL, ixCL, ixCCL, igD, nxD, normconv, wnamei, wC, wCC)
Definition: convert.t:1852
character(len=std_len) base_filename
Base file name for simulation output, which will be followed by a number.
Module with basic grid data structures.
Definition: mod_forest.t:2
integer itag
MPI recv send variables for AMR.
subroutine calc_grid(qunit, igrid, xC, xCC, xC_TMP, xCC_TMP, wC_TMP, wCC_TMP, normconv, ixCL, ixCCL, first)
Compute both corner as well as cell-centered values for output.
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl...
Definition: mod_physics.t:22
integer npe
The number of MPI tasks.
subroutine save_conntec(qunit, igrid, igonlevel)
Definition: convert.t:566
logical b0field
split magnetic field as background B0 field
logical, dimension(:), allocatable writelevel
subroutine unstructuredvtk(qunit)
Definition: convert.t:641
integer function nodenumbertec1d(i1, nx1, ig, igrid)
Definition: convert.t:614
Pointer to a tree_node.
Definition: mod_forest.t:6
integer, dimension(:), allocatable ng
number of grid blocks in domain per dimension, in array over levels
subroutine calc_x(igrid, xC, xCC)
computes cell corner (xC) and cell center (xCC) coordinates
procedure(aux_output), pointer usr_aux_output
Module with all the methods that users can customize in AMRVAC.
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 type_block_io
MPI type for IO: block excluding ghost cells.
integer, parameter unitconvert
subroutine punstructuredvtkb_mpi(qunit)
Definition: convert.t:2471
logical nocartesian
IO switches for conversion.
integer ierrmpi
A global MPI error return code.
integer ixm
the mesh range (within a block with ghost cells)
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
subroutine unstructuredvtk_mpi(qunit)
Definition: convert.t:1554
subroutine tecplot(qunit)
Definition: convert.t:364
logical image
output image
subroutine write_vtk(qunit, ixIL, ixCL, ixCCL, igrid, nc, np, nxD, nxCD, normconv, wnamei, xC, xCC, wC, wCC)
Definition: convert.t:1743
double precision, dimension(:), allocatable w_convert_factor
Conversion factors the primitive variables.
subroutine punstructuredvtk_mpi(qunit)
Definition: convert.t:1460
integer, dimension(:), allocatable, save morton_start
First Morton number per processor.
Definition: mod_forest.t:62
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
subroutine unstructuredvtkb(qunit)
Definition: convert.t:821
integer, parameter unitterm
Unit for standard output.
integer, dimension(:), allocatable, save morton_stop
Last Morton number per processor.
Definition: mod_forest.t:65
integer mype
The rank of the current MPI task.
double precision global_time
The global simulation time.
subroutine resettree_convert
Force the tree to desired level(s) from level_io(_min/_max) used for conversion to vtk output...
Definition: amrgrid.t:67
subroutine get_euv_spectra(qunit)
subroutine save_connvtk(qunit, igrid)
Definition: convert.t:1251
subroutine tecplot_mpi(qunit)
Definition: convert.t:1987
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
integer function nodenumbertec3d(i1, i2, i3, nx1, nx2, nx3, ig, igrid)
Definition: convert.t:632
double precision length_convert_factor
double precision, dimension(ndim) dxlevel
integer snapshotini
Resume from the snapshot with this index.
integer, parameter ndim
Number of spatial dimensions for grid variables.
subroutine generate_plotfile
Definition: convert.t:3
integer nwauxio
Number of auxiliary variables that are only included in output.
integer icomm
The MPI communicator.
integer, parameter nlevelshi
The maximum number of levels in the grid refinement.
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:47
integer, dimension(:,:), allocatable node
logical phys_energy
Solve energy equation or not.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer type_block_xc_io
MPI type for IO: cell corner (xc) or cell center (xcc) coordinates.
logical spectrum
output spectrum
logical autoconvert
If true, already convert to output format during the run.
integer function nodenumbertec2d(i1, i2, nx1, nx2, ig, igrid)
Definition: convert.t:623
procedure(special_convert), pointer usr_special_convert
Handles computations for coordinates and variables in output.