MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_convert_files.t
Go to the documentation of this file.
2 use mod_comm_lib, only: mpistop
3
4 implicit none
5 public
6
7contains
8
14 use mod_convert, only: convert_all
16
17 character(len=std_len) :: convert_type_elem
18 integer :: i
19
20 select case(convert_type)
21 case('tecplot','tecplotCC','tecline')
23 case('tecplotmpi','tecplotCCmpi','teclinempi')
25 case('vtu','vtuCC')
27 case('vtumpi','vtuCCmpi')
29 case('vtuB','vtuBCC','vtuBmpi','vtuBCCmpi')
31 case('vtuB64','vtuBCC64','vtuBmpi64','vtuBCCmpi64')
33 {^iftwod
34 case('vtuB23','vtuBCC23')
36 case('vtuBsym23','vtuBCCsym23')
38 }
39 case('pvtumpi','pvtuCCmpi')
41 case('pvtuBmpi','pvtuBCCmpi')
43 case('vtimpi','vtiCCmpi')
45 case('onegrid','onegridmpi')
47 case('oneblock','oneblockB')
49 case('EIvtiCCmpi','ESvtiCCmpi','SIvtiCCmpi','WIvtiCCmpi','EIvtuCCmpi','ESvtuCCmpi','SIvtuCCmpi','WIvtuCCmpi')
50 ! output synthetic euv emission
51 if (ndim==3 .and. associated(phys_te_images)) then
53 endif
54 case('dat_generic_mpi')
55 call convert_all()
56 case('user','usermpi')
57 if (.not. associated(usr_special_convert)) then
58 call mpistop("usr_special_convert not defined")
59 else
61 end if
62 case default
63 call mpistop("Error in generate_plotfile: Unknown convert_type")
64 end select
65
66 end subroutine generate_plotfile
67
68 subroutine oneblock(qunit)
69 ! this is for turning an AMR run into a single block
70 ! the data will be all on selected level level_io
71 ! this version should work for any dimension
72 ! only writes w_write selected 1:nw variables, also nwauxio
73 ! may use saveprim to switch to primitives
74 ! this version can not work on multiple CPUs
75 ! does not renormalize variables
76 ! header info differs from onegrid below
77 ! ASCII or binary output
78
82 use mod_physics
85 integer, intent(in) :: qunit
86
87 double precision:: normconv(0:nw+nwauxio)
88 integer :: Morton_no,igrid,ix^D,ig^D,level
89 integer, pointer :: ig_to_igrid(:^D&,:)
90 integer :: filenr,ncells^D,ncellx^D,jg^D,jig^D
91 integer :: iw,iiw,writenw,iwrite(1:nw+nwauxio),iigrid,idim
92 logical :: fileopen,writeblk(max_blocks)
93 logical :: patchw(ixG^T)
94 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
95 character(len=1024) :: outfilehead
96 character(len=80) :: filename
97
98 if(level_io<1)then
99 call mpistop('please specify level_io>0 for usage with oneblock')
100 end if
101
102 if(autoconvert)then
103 call mpistop('Set autoconvert=F and convert oneblock data manually')
104 end if
105
106 if(npe>1)then
107 if(mype==0) print *,'ONEBLOCK as yet to be parallelized'
108 call mpistop('npe>1, oneblock')
109 end if
110
111 ! only variables selected by w_write will be written out
112 normconv(0:nw+nwauxio)=one
113 normconv(0) = length_convert_factor
114 normconv(1:nw) = w_convert_factor
115 writenw=count(w_write(1:nw))+nwauxio
116 iiw=0
117 do iw =1,nw
118 if (.not.w_write(iw))cycle
119 iiw=iiw+1
120 iwrite(iiw)=iw
121 end do
122 if(nwauxio>0)then
123 do iw =nw+1,nw+nwauxio
124 iiw=iiw+1
125 iwrite(iiw)=iw
126 end do
127 end if
128
129 allocate(ig_to_igrid(ng^d(level_io),0:npe-1))
130 ig_to_igrid=-1
131 writeblk=.false.
132 do morton_no=morton_start(mype),morton_stop(mype)
133 igrid=sfc_to_igrid(morton_no)
134 level=node(plevel_,igrid)
135 ig^d=igrid_to_node(igrid,mype)%node%ig^d;
136 ig_to_igrid(ig^d,mype)=igrid
137 if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
138 *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
139 <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
140 writeblk(igrid)=.true.
141 end if
142 end do
143
144 call getheadernames(wnamei,xandwnamei,outfilehead)
145 ncells^d=0;
146 ncellx^d=ixmhi^d-ixmlo^d+1\
147 {do ig^d=1,ng^d(level_io)\}
148 igrid=ig_to_igrid(ig^d,mype)
149 if(writeblk(igrid)) go to 20
150 {end do\}
151 20 continue
152 jg^d=ig^d;
153 {
154 jig^dd=jg^dd;
155 do ig^d=1,ng^d(level_io)
156 jig^d=ig^d
157 igrid=ig_to_igrid(jig^dd,mype)
158 if(writeblk(igrid)) ncells^d=ncells^d+ncellx^d
159 end do
160 \}
161
162 do iigrid=1,igridstail; igrid=igrids(iigrid)
163 if(.not.writeblk(igrid)) cycle
164 block=>ps(igrid)
165 call alloc_state_output(igrid,ps1(igrid),ixg^ll)
166 ps1(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
167 if(nwauxio > 0) then
168 if (.not. associated(usr_aux_output)) then
169 call mpistop("usr_aux_output not defined")
170 else
171 call usr_aux_output(ixg^ll,ixm^ll^ladd1,ps1(igrid)%w,ps(igrid)%x,normconv)
172 end if
173 end if
174 if(saveprim) then
175 call phys_to_primitive(ixg^ll,ixg^ll^lsub2,ps1(igrid)%w,ps(igrid)%x)
176 if(b0field) ps1(igrid)%w(ixg^t,iw_mag(:))=ps1(igrid)%w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
177 else
178 if(b0field) then
179 ! add background magnetic field B0 to B
180 if(phys_total_energy) then
181 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) &
182 + sum(ps1(igrid)%w(ixg^t,iw_mag(:))*ps(igrid)%B0(ixg^t,:,0),dim=ndim+1)
183 end if
184 ps1(igrid)%w(ixg^t,iw_mag(:))=ps1(igrid)%w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
185 end if
186 end if
187 end do
188
189 master_cpu_open : if (mype == 0) then
190 inquire(qunit,opened=fileopen)
191 if (.not.fileopen) then
192 ! generate filename
193 filenr=snapshotini
194 if (autoconvert) filenr=snapshotnext
195 write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".blk"
196 select case(convert_type)
197 case("oneblock")
198 open(qunit,file=filename,status='unknown')
199 write(qunit,*) trim(outfilehead)
200 write(qunit,*) ncells^d
201 write(qunit,*) real(global_time*time_convert_factor)
202 case("oneblockB")
203 open(qunit,file=filename,form='unformatted',status='unknown')
204 write(qunit) outfilehead
205 write(qunit) ncells^d
206 write(qunit) real(global_time*time_convert_factor)
207 end select
208 end if
209 end if master_cpu_open
210
211 {^ifthreed
212 do ig3=1,ng3(level_io)
213 do ix3=ixmlo3,ixmhi3}
214
215 {^nooned
216 do ig2=1,ng2(level_io)
217 do ix2=ixmlo2,ixmhi2}
218
219 do ig1=1,ng1(level_io)
220 igrid=ig_to_igrid(ig^d,mype)
221 if(.not.writeblk(igrid)) cycle
222 do ix1=ixmlo1,ixmhi1
223 master_write : if(mype==0) then
224 select case(convert_type)
225 case("oneblock")
226 write(qunit,fmt="(100(e14.6))") &
227 ps(igrid)%x(ix^d,1:ndim)*normconv(0),&
228 (ps1(igrid)%w(ix^d,iwrite(iw))*normconv(iwrite(iw)),iw=1,writenw)
229 case("oneblockB")
230 write(qunit) real(ps(igrid)%x(ix^d,1:ndim)*normconv(0)),&
231 (real(ps1(igrid)%w(ix^d,iwrite(iw))*normconv(iwrite(iw))),iw=1,writenw)
232 end select
233 end if master_write
234 end do
235 end do
236 {^nooned
237 end do
238 end do}
239 {^ifthreed
240 end do
241 end do}
242
243 close(qunit)
244
245 end subroutine oneblock
246
247 subroutine onegrid(qunit)
248 ! this is for turning an AMR run into a single grid
249 ! this version should work for any dimension, can be in parallel
250 ! in 1D, should behave much like oneblock, except for header info
251
252 ! only writes all 1:nw variables, no nwauxio
253 ! may use saveprim to switch to primitives
254 ! this version can work on multiple CPUs
255 ! does not renormalize variables
256 ! ASCII output
257
260 use mod_physics
262 integer, intent(in) :: qunit
263
264 double precision :: w_recv(ixG^T,1:nw),x_recv(ixG^T,1:ndim)
265 integer :: itag,Morton_no,igrid,ix^D,iw
266 integer :: filenr
267 !.. MPI variables ..
268 integer :: igrid_recv,ipe
269 integer, allocatable :: intstatus(:,:)
270 logical :: fileopen
271 character(len=80) :: filename
272 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
273 character(len=1024) :: outfilehead
274
275
276 if(nwauxio>0)then
277 if(mype==0) print *,'ONEGRID to be used without nwauxio'
278 call mpistop('nwauxio>0, onegrid')
279 end if
280
281 if(saveprim)then
282 if(mype==0.and.nwaux>0) print *,'warning: ONEGRID used with saveprim, check auxiliaries'
283 end if
284
285 master_cpu_open : if (mype == 0) then
286 call getheadernames(wnamei,xandwnamei,outfilehead)
287 write(outfilehead,'(a)') "#"//" "//trim(outfilehead)
288 inquire(qunit,opened=fileopen)
289 if (.not.fileopen) then
290 ! generate filename
291 filenr=snapshotini
292 if (autoconvert) filenr=snapshotnext
293 write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".blk"
294 open(qunit,file=filename,status='unknown')
295 end if
296 write(qunit,"(a)")outfilehead
297 write(qunit,"(i7)") ( {^d&(ixmhi^d-ixmlo^d+1)*} )*(morton_stop(npe-1)-morton_start(0)+1)
298 end if master_cpu_open
299
300 do morton_no=morton_start(mype),morton_stop(mype)
301 igrid=sfc_to_igrid(morton_no)
302 if(saveprim) call phys_to_primitive(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
303 if(mype/=0)then
304 itag=morton_no
305 call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
306 call mpi_send(ps(igrid)%x,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
307 itag=igrid
308 call mpi_send(ps(igrid)%w,1,type_block_io, 0,itag,icomm,ierrmpi)
309 else
310 {do ix^db=ixmlo^db,ixmhi^db\}
311 do iw=1,nw
312 if( dabs(ps(igrid)%w(ix^d,iw)) < 1.0d-32 ) ps(igrid)%w(ix^d,iw) = zero
313 end do
314 write(qunit,fmt="(100(e14.6))") ps(igrid)%x(ix^d,1:ndim),ps(igrid)%w(ix^d,1:nw)
315 {end do\}
316 end if
317 end do
318
319 if(mype==0.and.npe>1) allocate(intstatus(mpi_status_size,1))
320
321 manycpu : if (npe>1) then
322 if (mype==0) then
323 loop_cpu : do ipe =1, npe-1
324 loop_morton : do morton_no=morton_start(ipe),morton_stop(ipe)
325 itag=morton_no
326 call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
327 call mpi_recv(x_recv,1,type_block_xcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
328 itag=igrid_recv
329 call mpi_recv(w_recv,1,type_block_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
330 {do ix^db=ixmlo^db,ixmhi^db\}
331 do iw=1,nw
332 if( dabs(ps(igrid)%w(ix^d,iw)) < smalldouble ) ps(igrid)%w(ix^d,iw) = zero
333 end do
334 write(qunit,fmt="(100(e14.6))") x_recv(ix^d,1:ndim),w_recv(ix^d,1:nw)
335 {end do\}
336 end do loop_morton
337 end do loop_cpu
338 end if
339 end if manycpu
340
341 if (npe>1) then
342 call mpi_barrier(icomm,ierrmpi)
343 if(mype==0)deallocate(intstatus)
344 end if
345
346 if(mype==0) close(qunit)
347 end subroutine onegrid
348
349 subroutine tecplot(qunit)
350
351 ! output for tecplot (ASCII format)
352 ! not parallel, uses calc_grid to compute nwauxio variables
353 ! allows renormalizing using convert factors
354
357
358 integer, intent(in) :: qunit
359
360 double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
361 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
362 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
363 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
364 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
365 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
366 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
367 double precision, dimension(0:nw+nwauxio) :: normconv
368 integer:: igrid,iigrid,level,igonlevel,iw,idim,ix^D
369 integer:: NumGridsOnLevel(1:nlevelshi)
370 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,ixC^L,ixCC^L
371 integer :: nodes, elems
372 integer :: filenr
373 logical :: fileopen,first
374 character(len=80) :: filename
375 !!! possible length conflict
376 character(len=1024) :: tecplothead
377 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
378 character(len=1024) :: outfilehead
379
380 if(npe>1)then
381 if(mype==0) print *,'tecplot not parallel, use tecplotmpi'
382 call mpistop('npe>1, tecplot')
383 end if
384
385 if(nw/=count(w_write(1:nw)))then
386 if(mype==0) print *,'tecplot does not use w_write=F'
387 call mpistop('w_write, tecplot')
388 end if
389
390 if(nocartesian)then
391 if(mype==0) print *,'tecplot with nocartesian'
392 end if
393
394 inquire(qunit,opened=fileopen)
395 if(.not.fileopen) then
396 ! generate filename
397 filenr=snapshotini
398 if (autoconvert) filenr=snapshotnext
399 write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".plt"
400 open(qunit,file=filename,status='unknown')
401 end if
402
403 call getheadernames(wnamei,xandwnamei,outfilehead)
404
405 write(tecplothead,'(a)') "VARIABLES = "//trim(outfilehead)
406 write(qunit,'(a)') tecplothead(1:len_trim(tecplothead))
407
408 numgridsonlevel(1:nlevelshi)=0
409 do level=levmin,levmax
410 numgridsonlevel(level)=0
411 do iigrid=1,igridstail; igrid=igrids(iigrid);
412 if (node(plevel_,igrid)/=level) cycle
413 numgridsonlevel(level)=numgridsonlevel(level)+1
414 end do
415 end do
416
417 nx^d=ixmhi^d-ixmlo^d+1;
418 nxc^d=nx^d+1;
419
420 {^ifoned
421 if(convert_type=='tecline') then
422 nodes=0
423 elems=0
424 do level=levmin,levmax
425 nodes=nodes + numgridsonlevel(level)*{nxc^d*}
426 elems=elems + numgridsonlevel(level)*{nx^d*}
427 end do
428
429 write(qunit,"(a,i7,a,1pe12.5,a)") &
430 'ZONE T="all levels", I=',elems, &
431 ', SOLUTIONTIME=',global_time*time_convert_factor,', F=POINT'
432
433 igonlevel=0
434 do iigrid=1,igridstail; igrid=igrids(iigrid);
435 block=>ps(igrid)
436 call calc_x(igrid,xc,xcc)
437 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,ixc^l,ixcc^l,.true.)
438 {do ix^db=ixccmin^db,ixccmax^db\}
439 x_tec(1:ndim)=xcc_tmp(ix^d,1:ndim)*normconv(0)
440 w_tec(1:nw+nwauxio)=wcc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
441 write(qunit,fmt="(100(e24.16))") x_tec, w_tec
442 {end do\}
443 end do
444 close(qunit)
445 else
446 }
447 do level=levmin,levmax
448 nodesonlevel=numgridsonlevel(level)*{nxc^d*}
449 elemsonlevel=numgridsonlevel(level)*{nx^d*}
450 ! for all tecplot variants coded up here, we let the TECPLOT ZONES coincide
451 ! with the AMR grid LEVEL. Other options would be
452 ! let each grid define a zone: inefficient for TECPLOT internal workings
453 ! hence not implemented
454 ! let entire octree define 1 zone: no difference in interpolation
455 ! properties across TECPLOT zones detected as yet, hence not done
456 select case(convert_type)
457 case('tecplot')
458 ! in this option, we store the corner coordinates, as well as the corner
459 ! values of all variables (obtained by averaging). This allows POINT packaging,
460 ! and thus we can save full grid info by using one call to calc_grid
461 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
462 'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
463 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
464 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
465 do iigrid=1,igridstail; igrid=igrids(iigrid);
466 if (node(plevel_,igrid)/=level) cycle
467 block=>ps(igrid)
468 call calc_x(igrid,xc,xcc)
469 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
470 ixc^l,ixcc^l,.true.)
471 {do ix^db=ixcmin^db,ixcmax^db\}
472 x_tec(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0)
473 w_tec(1:nw+nwauxio)=wc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
474 write(qunit,fmt="(100(e14.6))") x_tec, w_tec
475 {end do\}
476 end do
477 case('tecplotCC')
478 ! in this option, we store the corner coordinates, and the cell center
479 ! values of all variables. Due to this mix of corner/cell center, we must
480 ! use BLOCK packaging, and thus we have enormous overhead by using
481 ! calc_grid repeatedly to merely fill values of cell corner coordinates
482 ! and cell center values per dimension, per variable
483 if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
484 if(nw+nwauxio==1)then
485 ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
486 ! and just set [ndim+1]
487 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
488 'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
489 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
490 ndim+1,']=CELLCENTERED), ZONETYPE=', &
491 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
492 else
493 if(ndim+nw+nwauxio<10) then
494 ! difference only in length of integer format specification for ndim+nw+nwauxio
495 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
496 'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
497 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
498 ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
499 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
500 else
501 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
502 'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
503 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
504 ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
505 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
506 end if
507 end if
508 do idim=1,ndim
509 first=(idim==1)
510 do iigrid=1,igridstail; igrid=igrids(iigrid);
511 if (node(plevel_,igrid)/=level) cycle
512 block=>ps(igrid)
513 call calc_x(igrid,xc,xcc)
514 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
515 ixc^l,ixcc^l,first)
516 write(qunit,fmt="(100(e14.6))") xc_tmp(ixc^s,idim)*normconv(0)
517 end do
518 end do
519 do iw=1,nw+nwauxio
520 do iigrid=1,igridstail; igrid=igrids(iigrid);
521 if (node(plevel_,igrid)/=level) cycle
522 block=>ps(igrid)
523 call calc_x(igrid,xc,xcc)
524 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
525 ixc^l,ixcc^l,.true.)
526 write(qunit,fmt="(100(e14.6))") wcc_tmp(ixcc^s,iw)*normconv(iw)
527 enddo
528 enddo
529 case default
530 call mpistop('no such tecplot type')
531 end select
532 igonlevel=0
533 do iigrid=1,igridstail; igrid=igrids(iigrid);
534 if (node(plevel_,igrid)/=level) cycle
535 block=>ps(igrid)
536 igonlevel=igonlevel+1
537 call save_conntec(qunit,igrid,igonlevel)
538 end do
539 end do
540 {^ifoned endif}
541
542 close(qunit)
543
544 end subroutine tecplot
545
546 subroutine save_conntec(qunit,igrid,igonlevel)
547
548 ! this saves the basic line, quad and brick connectivity,
549 ! as used by TECPLOT file outputs for unstructured grid
551
552 integer, intent(in) :: qunit, igrid, igonlevel
553
554 integer :: nx^D, nxC^D, ix^D
555
556 nx^d=ixmhi^d-ixmlo^d+1;
557 nxc^d=nx^d+1;
558
559 ! connectivity list
560 {do ix^db=1,nx^db\}
561 {^ifthreed
562 ! basic brick connectivity
563 write(qunit,'(8(i7,1x))') &
564 nodenumbertec3d(ix1, ix2-1,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
565 nodenumbertec3d(ix1+1,ix2-1,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
566 nodenumbertec3d(ix1+1,ix2 ,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
567 nodenumbertec3d(ix1 ,ix2 ,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
568 nodenumbertec3d(ix1 ,ix2-1,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
569 nodenumbertec3d(ix1+1,ix2-1,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
570 nodenumbertec3d(ix1+1,ix2 ,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
571 nodenumbertec3d(ix1 ,ix2 ,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid)
572 }
573 {^iftwod
574 ! basic quadrilateral connectivity
575 write(qunit,'(4(i7,1x))') &
576 nodenumbertec2d(ix1, ix2-1,nxc1,nxc2,igonlevel,igrid),&
577 nodenumbertec2d(ix1+1,ix2-1,nxc1,nxc2,igonlevel,igrid),&
578 nodenumbertec2d(ix1+1,ix2 ,nxc1,nxc2,igonlevel,igrid),&
579 nodenumbertec2d(ix1 ,ix2 ,nxc1,nxc2,igonlevel,igrid)
580 }
581 {^ifoned
582 ! basic line connectivity
583 write(qunit,'(2(i7,1x))') nodenumbertec1d(ix1,nxc1,igonlevel,igrid),&
584 nodenumbertec1d(ix1+1,nxc1,igonlevel,igrid)
585 }
586 {end do\}
587
588 end subroutine save_conntec
589
590 integer function nodenumbertec1d(i1,nx1,ig,igrid)
591 use mod_comm_lib, only: mpistop
592 integer, intent(in):: i1,nx1,ig,igrid
593
594 nodenumbertec1d=i1+(ig-1)*nx1
595 if(nodenumbertec1d>9999999)call mpistop("too large nodenumber")
596
597 end function nodenumbertec1d
598
599 integer function nodenumbertec2d(i1,i2,nx1,nx2,ig,igrid)
600
601 integer, intent(in):: i1,i2,nx1,nx2,ig,igrid
602
603 nodenumbertec2d=i1+i2*nx1+(ig-1)*nx1*nx2
604 if(nodenumbertec2d>9999999)call mpistop("too large nodenumber")
605
606 end function nodenumbertec2d
607
608 integer function nodenumbertec3d(i1,i2,i3,nx1,nx2,nx3,ig,igrid)
609
610 integer, intent(in):: i1,i2,i3,nx1,nx2,nx3,ig,igrid
611
612 nodenumbertec3d=i1+i2*nx1+i3*nx1*nx2+(ig-1)*nx1*nx2*nx3
613 if(nodenumbertec3d>9999999)call mpistop("too large nodenumber")
614
615 end function nodenumbertec3d
616
617 subroutine unstructuredvtk(qunit)
618
619 ! output for vtu format to paraview
620 ! not parallel, uses calc_grid to compute nwauxio variables
621 ! allows renormalizing using convert factors
624
625 integer, intent(in) :: qunit
626
627 double precision :: x_VTK(1:3)
628 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
629 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
630 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
631 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
632 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
633 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
634 double precision, dimension(0:nw+nwauxio) :: normconv
635 integer:: igrid,iigrid,level,igonlevel,icel,ixC^L,ixCC^L,iw
636 integer:: NumGridsOnLevel(1:nlevelshi)
637 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,VTK_type,ix^D
638 integer :: filenr
639 character(len=80):: filename
640 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
641 character(len=1024) :: outfilehead
642 logical :: fileopen
643
644 if(npe>1)then
645 if(mype==0) print *,'unstructuredvtk not parallel, use vtumpi'
646 call mpistop('npe>1, unstructuredvtk')
647 end if
648
649 inquire(qunit,opened=fileopen)
650 if(.not.fileopen)then
651 ! generate filename
652 filenr=snapshotini
653 if (autoconvert) filenr=snapshotnext
654 write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
655 ! Open the file for the header part
656 open(qunit,file=filename,status='unknown')
657 end if
658
659 call getheadernames(wnamei,xandwnamei,outfilehead)
660
661 ! generate xml header
662 write(qunit,'(a)')'<?xml version="1.0"?>'
663 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
664 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
665 write(qunit,'(a)')'<UnstructuredGrid>'
666 write(qunit,'(a)')'<FieldData>'
667 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
668 'NumberOfTuples="1" format="ascii">'
669 write(qunit,*) dble(dble(global_time)*time_convert_factor)
670 write(qunit,'(a)')'</DataArray>'
671 write(qunit,'(a)')'</FieldData>'
672
673 ! number of cells, number of corner points, per grid.
674 nx^d=ixmhi^d-ixmlo^d+1;
675 nxc^d=nx^d+1;
676 nc={nx^d*}
677 np={nxc^d*}
678
679 ! Note: using the w_write, writelevel, writespshift
680 ! we can clip parts of the grid away, select variables, levels etc.
681 do level=levmin,levmax
682 if (writelevel(level)) then
683 do iigrid=1,igridstail; igrid=igrids(iigrid);
684 if (node(plevel_,igrid)/=level) cycle
685 block=>ps(igrid)
686 ! only output a grid when fully within clipped region selected
687 ! by writespshift array
688 if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
689 *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
690 <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
691 call calc_x(igrid,xc,xcc)
692 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
693 ixc^l,ixcc^l,.true.)
694 select case(convert_type)
695 case('vtu')
696 ! we write out every grid as one VTK PIECE
697 write(qunit,'(a,i7,a,i7,a)') &
698 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
699 write(qunit,'(a)')'<PointData>'
700 do iw=1,nw+nwauxio
701 if(iw<=nw) then
702 if(.not.w_write(iw)) cycle
703 end if
704 write(qunit,'(a,a,a)')&
705 '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
706 write(qunit,'(200(1pe14.6))') {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
707 write(qunit,'(a)')'</DataArray>'
708 end do
709 write(qunit,'(a)')'</PointData>'
710 write(qunit,'(a)')'<Points>'
711 write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
712 ! write cell corner coordinates in a backward dimensional loop, always 3D output
713 {do ix^db=ixcmin^db,ixcmax^db \}
714 x_vtk(1:3)=zero;
715 x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
716 write(qunit,'(3(1pe14.6))') x_vtk
717 {end do \}
718 write(qunit,'(a)')'</DataArray>'
719 write(qunit,'(a)')'</Points>'
720 case('vtuCC')
721 ! we write out every grid as one VTK PIECE
722 write(qunit,'(a,i7,a,i7,a)') &
723 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
724 write(qunit,'(a)')'<CellData>'
725 do iw=1,nw+nwauxio
726 if(iw<=nw) then
727 if(.not.w_write(iw)) cycle
728 end if
729 write(qunit,'(a,a,a)')&
730 '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
731 write(qunit,'(200(1pe14.6))') {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
732 write(qunit,'(a)')'</DataArray>'
733 end do
734 write(qunit,'(a)')'</CellData>'
735 write(qunit,'(a)')'<Points>'
736 write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
737 ! write cell corner coordinates in a backward dimensional loop, always 3D output
738 {do ix^db=ixcmin^db,ixcmax^db \}
739 x_vtk(1:3)=zero;
740 x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
741 write(qunit,'(3(1pe14.6))') x_vtk
742 {end do \}
743 write(qunit,'(a)')'</DataArray>'
744 write(qunit,'(a)')'</Points>'
745 end select
746
747 write(qunit,'(a)')'<Cells>'
748 ! connectivity part
749 write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
750 call save_connvtk(qunit,igrid)
751 write(qunit,'(a)')'</DataArray>'
752
753 ! offsets data array
754 write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
755 do icel=1,nc
756 write(qunit,'(i7)') icel*(2**^nd)
757 end do
758 write(qunit,'(a)')'</DataArray>'
759
760 ! VTK cell type data array
761 write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
762 ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
763 {^ifoned vtk_type=3 \}
764 {^iftwod vtk_type=8 \}
765 {^ifthreed vtk_type=11 \}
766 do icel=1,nc
767 write(qunit,'(i2)') vtk_type
768 end do
769 write(qunit,'(a)')'</DataArray>'
770
771 write(qunit,'(a)')'</Cells>'
772
773 write(qunit,'(a)')'</Piece>'
774 end if
775 end do
776 end if
777 end do
778
779 write(qunit,'(a)')'</UnstructuredGrid>'
780 write(qunit,'(a)')'</VTKFile>'
781 close(qunit)
782
783 end subroutine unstructuredvtk
784
785 subroutine unstructuredvtkb(qunit)
786
787 ! output for vtu format to paraview, binary version output
788 ! not parallel, uses calc_grid to compute nwauxio variables
789 ! allows renormalizing using convert factors
793
794 integer, intent(in) :: qunit
795
796 double precision :: x_VTK(1:3)
797 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
798 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
799 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
800 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
801 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio):: wC_TMP
802 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
803 double precision :: normconv(0:nw+nwauxio)
804 integer, allocatable :: intstatus(:,:)
805 integer*8 :: offset
806 integer :: itag,ipe,igrid,level,icel,ixC^L,ixCC^L,Morton_no,Morton_length
807 integer :: nx^D,nxC^D,nc,np,VTK_type,ix^D,filenr
808 integer:: k,iw
809 integer:: length,lengthcc,length_coords,length_conn,length_offsets
810 character:: buf
811 character(len=80):: filename
812 character(len=19):: offset_char
813 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
814 character(len=1024) :: outfilehead
815 logical :: fileopen,cell_corner=.false.
816 logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
817
818 normconv=one
819 morton_length=morton_stop(npe-1)-morton_start(0)+1
820 allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
821 allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
822 morton_aim=.false.
823 morton_aim_p=.false.
824 do morton_no=morton_start(mype),morton_stop(mype)
825 igrid=sfc_to_igrid(morton_no)
826 level=node(plevel_,igrid)
827 ! we can clip parts of the grid away, select variables, levels etc.
828 if(writelevel(level)) then
829 ! only output a grid when fully within clipped region selected
830 ! by writespshift array
831 if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
832 *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
833 <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
834 morton_aim_p(morton_no)=.true.
835 end if
836 end if
837 end do
838 call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
840 select case(convert_type)
841 case('vtuB','vtuBmpi')
842 cell_corner=.true.
843 case('vtuBCC','vtuBCCmpi')
844 cell_corner=.false.
845 end select
846 if (mype /= 0) then
847 do morton_no=morton_start(mype),morton_stop(mype)
848 if(.not. morton_aim(morton_no)) cycle
849 igrid=sfc_to_igrid(morton_no)
850 call calc_x(igrid,xc,xcc)
851 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
852 ixc^l,ixcc^l,.true.)
853 itag=morton_no
854 call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
855 if(cell_corner) then
856 call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
857 else
858 call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
859 end if
860 end do
861
862 else
863 ! mype==0
864 offset=0
865 inquire(qunit,opened=fileopen)
866 if(.not.fileopen)then
867 ! generate filename
868 filenr=snapshotini
869 if (autoconvert) filenr=snapshotnext
870 write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
871 ! Open the file for the header part
872 open(qunit,file=filename,status='replace')
873 end if
874 call getheadernames(wnamei,xandwnamei,outfilehead)
875 ! generate xml header
876 write(qunit,'(a)')'<?xml version="1.0"?>'
877 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
878 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
879 write(qunit,'(a)')'<UnstructuredGrid>'
880 write(qunit,'(a)')'<FieldData>'
881 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
882 'NumberOfTuples="1" format="ascii">'
883 write(qunit,*) real(global_time*time_convert_factor)
884 write(qunit,'(a)')'</DataArray>'
885 write(qunit,'(a)')'</FieldData>'
886
887 ! number of cells, number of corner points, per grid.
888 nx^d=ixmhi^d-ixmlo^d+1;
889 nxc^d=nx^d+1;
890 nc={nx^d*}
891 np={nxc^d*}
892 length=np*size_real
893 lengthcc=nc*size_real
894 length_coords=3*length
895 length_conn=2**^nd*size_int*nc
896 length_offsets=nc*size_int
897
898 ! Note: using the w_write, writelevel, writespshift
899 do morton_no=morton_start(0),morton_stop(0)
900 if(.not. morton_aim(morton_no)) cycle
901 if(cell_corner) then
902 ! we write out every grid as one VTK PIECE
903 write(qunit,'(a,i7,a,i7,a)') &
904 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
905 write(qunit,'(a)')'<PointData>'
906 do iw=1,nw+nwauxio
907 if(iw<=nw) then
908 if(.not.w_write(iw)) cycle
909 endif
910 write(offset_char,'(i19)') offset
911 write(qunit,'(a,a,a,a,a)')&
912 '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
913 '" format="appended" offset="',trim(adjustl(offset_char)),'">'
914 write(qunit,'(a)')'</DataArray>'
915 offset=offset+length+size_int
916 end do
917 write(qunit,'(a)')'</PointData>'
918 write(qunit,'(a)')'<Points>'
919 write(offset_char,'(i19)') offset
920 write(qunit,'(a,a,a)') &
921 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
922 ! write cell corner coordinates in a backward dimensional loop, always 3D output
923 offset=offset+length_coords+size_int
924 write(qunit,'(a)')'</Points>'
925 else
926 ! we write out every grid as one VTK PIECE
927 write(qunit,'(a,i7,a,i7,a)') &
928 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
929 write(qunit,'(a)')'<CellData>'
930 do iw=1,nw+nwauxio
931 if(iw<=nw) then
932 if(.not.w_write(iw)) cycle
933 end if
934 write(offset_char,'(i19)') offset
935 write(qunit,'(a,a,a,a,a)')&
936 '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
937 '" format="appended" offset="',trim(adjustl(offset_char)),'">'
938 write(qunit,'(a)')'</DataArray>'
939 offset=offset+lengthcc+size_int
940 end do
941 write(qunit,'(a)')'</CellData>'
942 write(qunit,'(a)')'<Points>'
943 write(offset_char,'(i19)') offset
944 write(qunit,'(a,a,a)') &
945 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
946 ! write cell corner coordinates in a backward dimensional loop, always 3D output
947 offset=offset+length_coords+size_int
948 write(qunit,'(a)')'</Points>'
949 end if
950 write(qunit,'(a)')'<Cells>'
951 ! connectivity part
952 write(offset_char,'(i19)') offset
953 write(qunit,'(a,a,a)')&
954 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
955 offset=offset+length_conn+size_int
956 ! offsets data array
957 write(offset_char,'(i19)') offset
958 write(qunit,'(a,a,a)') &
959 '<DataArray type="Int32" Name="offsets" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
960 offset=offset+length_offsets+size_int
961 ! VTK cell type data array
962 write(offset_char,'(i19)') offset
963 write(qunit,'(a,a,a)') &
964 '<DataArray type="Int32" Name="types" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
965 offset=offset+size_int+nc*size_int
966 write(qunit,'(a)')'</Cells>'
967 write(qunit,'(a)')'</Piece>'
968 end do
969 ! write metadata communicated from other processors
970 if(npe>1)then
971 do ipe=1, npe-1
972 do morton_no=morton_start(ipe),morton_stop(ipe)
973 if(.not. morton_aim(morton_no)) cycle
974 if(cell_corner) then
975 ! we write out every grid as one VTK PIECE
976 write(qunit,'(a,i7,a,i7,a)') &
977 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
978 write(qunit,'(a)')'<PointData>'
979 do iw=1,nw+nwauxio
980 if(iw<=nw) then
981 if(.not.w_write(iw)) cycle
982 end if
983 write(offset_char,'(i19)') offset
984 write(qunit,'(a,a,a,a,a)')&
985 '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
986 '" format="appended" offset="',trim(adjustl(offset_char)),'">'
987 write(qunit,'(a)')'</DataArray>'
988 offset=offset+length+size_int
989 end do
990 write(qunit,'(a)')'</PointData>'
991 write(qunit,'(a)')'<Points>'
992 write(offset_char,'(i19)') offset
993 write(qunit,'(a,a,a)') &
994 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
995 ! write cell corner coordinates in a backward dimensional loop, always 3D output
996 offset=offset+length_coords+size_int
997 write(qunit,'(a)')'</Points>'
998 else
999 ! we write out every grid as one VTK PIECE
1000 write(qunit,'(a,i7,a,i7,a)') &
1001 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1002 write(qunit,'(a)')'<CellData>'
1003 do iw=1,nw+nwauxio
1004 if(iw<=nw) then
1005 if(.not.w_write(iw)) cycle
1006 end if
1007 write(offset_char,'(i19)') offset
1008 write(qunit,'(a,a,a,a,a)')&
1009 '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
1010 '" format="appended" offset="',trim(adjustl(offset_char)),'">'
1011 write(qunit,'(a)')'</DataArray>'
1012 offset=offset+lengthcc+size_int
1013 end do
1014 write(qunit,'(a)')'</CellData>'
1015 write(qunit,'(a)')'<Points>'
1016 write(offset_char,'(i19)') offset
1017 write(qunit,'(a,a,a)') &
1018 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
1019 ! write cell corner coordinates in a backward dimensional loop, always 3D output
1020 offset=offset+length_coords+size_int
1021 write(qunit,'(a)')'</Points>'
1022 end if
1023 write(qunit,'(a)')'<Cells>'
1024 ! connectivity part
1025 write(offset_char,'(i19)') offset
1026 write(qunit,'(a,a,a)')&
1027 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
1028 offset=offset+length_conn+size_int
1029 ! offsets data array
1030 write(offset_char,'(i19)') offset
1031 write(qunit,'(a,a,a)') &
1032 '<DataArray type="Int32" Name="offsets" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
1033 offset=offset+length_offsets+size_int
1034 ! VTK cell type data array
1035 write(offset_char,'(i19)') offset
1036 write(qunit,'(a,a,a)') &
1037 '<DataArray type="Int32" Name="types" format="appended" offset="',trim(adjustl(offset_char)),'"/>'
1038 offset=offset+size_int+nc*size_int
1039 write(qunit,'(a)')'</Cells>'
1040 write(qunit,'(a)')'</Piece>'
1041 end do
1042 end do
1043 end if
1044
1045 write(qunit,'(a)')'</UnstructuredGrid>'
1046 write(qunit,'(a)')'<AppendedData encoding="raw">'
1047 close(qunit)
1048 open(qunit,file=filename,access='stream',form='unformatted',position='append')
1049 buf='_'
1050 write(qunit) trim(buf)
1051
1052 do morton_no=morton_start(0),morton_stop(0)
1053 if(.not. morton_aim(morton_no)) cycle
1054 igrid=sfc_to_igrid(morton_no)
1055 call calc_x(igrid,xc,xcc)
1056 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1057 ixc^l,ixcc^l,.true.)
1058 do iw=1,nw+nwauxio
1059 if(iw<=nw) then
1060 if(.not.w_write(iw)) cycle
1061 end if
1062 if(cell_corner) then
1063 write(qunit) length
1064 write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
1065 else
1066 write(qunit) lengthcc
1067 write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
1068 end if
1069 end do
1070
1071 write(qunit) length_coords
1072 {do ix^db=ixcmin^db,ixcmax^db \}
1073 x_vtk(1:3)=zero;
1074 x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1075 do k=1,3
1076 write(qunit) real(x_vtk(k))
1077 end do
1078 {end do \}
1079
1080 write(qunit) length_conn
1081 {do ix^db=1,nx^db\}
1082 {^ifoned write(qunit)ix1-1,ix1 \}
1083 {^iftwod
1084 write(qunit)(ix2-1)*nxc1+ix1-1, &
1085 (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1086 \}
1087 {^ifthreed
1088 write(qunit)&
1089 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1090 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1091 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1092 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1093 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1094 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1095 ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1096 ix3*nxc2*nxc1+ ix2*nxc1+ix1
1097 \}
1098 {end do\}
1099
1100 write(qunit) length_offsets
1101 do icel=1,nc
1102 write(qunit) icel*(2**^nd)
1103 end do
1104
1105 {^ifoned vtk_type=3 \}
1106 {^iftwod vtk_type=8 \}
1107 {^ifthreed vtk_type=11 \}
1108 write(qunit) size_int*nc
1109 do icel=1,nc
1110 write(qunit) vtk_type
1111 end do
1112 end do
1113 allocate(intstatus(mpi_status_size,1))
1114 if(npe>1)then
1115 ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d;
1116 ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
1117 do ipe=1, npe-1
1118 do morton_no=morton_start(ipe),morton_stop(ipe)
1119 if(.not. morton_aim(morton_no)) cycle
1120 itag=morton_no
1121 call mpi_recv(xc_tmp,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1122 if(cell_corner) then
1123 call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1124 else
1125 call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1126 end if
1127 do iw=1,nw+nwauxio
1128 if(iw<=nw) then
1129 if(.not.w_write(iw)) cycle
1130 end if
1131 if(cell_corner) then
1132 write(qunit) length
1133 write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
1134 else
1135 write(qunit) lengthcc
1136 write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
1137 end if
1138 end do
1139 write(qunit) length_coords
1140 {do ix^db=ixcmin^db,ixcmax^db \}
1141 x_vtk(1:3)=zero;
1142 x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1143 do k=1,3
1144 write(qunit) real(x_vtk(k))
1145 end do
1146 {end do \}
1147 write(qunit) length_conn
1148 {do ix^db=1,nx^db\}
1149 {^ifoned write(qunit)ix1-1,ix1 \}
1150 {^iftwod
1151 write(qunit)(ix2-1)*nxc1+ix1-1, &
1152 (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1153 \}
1154 {^ifthreed
1155 write(qunit)&
1156 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1157 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1158 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1159 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1160 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1161 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1162 ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1163 ix3*nxc2*nxc1+ ix2*nxc1+ix1
1164 \}
1165 {end do\}
1166 write(qunit) length_offsets
1167 do icel=1,nc
1168 write(qunit) icel*(2**^nd)
1169 end do
1170 {^ifoned vtk_type=3 \}
1171 {^iftwod vtk_type=8 \}
1172 {^ifthreed vtk_type=11 \}
1173 write(qunit) size_int*nc
1174 do icel=1,nc
1175 write(qunit) vtk_type
1176 end do
1177 end do
1178 end do
1179 end if
1180 close(qunit)
1181 open(qunit,file=filename,status='unknown',form='formatted',position='append')
1182 write(qunit,'(a)')'</AppendedData>'
1183 write(qunit,'(a)')'</VTKFile>'
1184 close(qunit)
1185 deallocate(intstatus)
1186 end if
1187
1188 deallocate(morton_aim,morton_aim_p)
1189 if (npe>1) then
1190 call mpi_barrier(icomm,ierrmpi)
1191 end if
1192
1193 end subroutine unstructuredvtkb
1194
1195 subroutine unstructuredvtkb64(qunit)
1196 ! output for vtu format to paraview, binary version output
1197 ! not parallel, uses calc_grid to compute nwauxio variables
1198 ! allows renormalizing using convert factors
1202
1203 integer, intent(in) :: qunit
1204
1205 double precision :: x_VTK(1:3)
1206 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
1207 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
1208 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1209 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1210 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio):: wC_TMP
1211 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
1212 double precision :: normconv(0:nw+nwauxio)
1213 integer, allocatable :: intstatus(:,:)
1214 integer*8 :: offset
1215 integer :: itag,ipe,igrid,level,icel,ixC^L,ixCC^L,Morton_no,Morton_length
1216 integer :: nx^D,nxC^D,nc,np,VTK_type,ix^D,filenr
1217 integer:: k,iw
1218 integer:: length,lengthcc,length_coords,length_conn,length_offsets
1219 character:: buf
1220 character(len=80):: filename
1221 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1222 character(len=1024) :: outfilehead
1223 logical :: fileopen,cell_corner=.false.
1224 logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
1225
1226 normconv=one
1227 morton_length=morton_stop(npe-1)-morton_start(0)+1
1228 allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
1229 allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
1230 morton_aim=.false.
1231 morton_aim_p=.false.
1232 do morton_no=morton_start(mype),morton_stop(mype)
1233 igrid=sfc_to_igrid(morton_no)
1234 level=node(plevel_,igrid)
1235 ! we can clip parts of the grid away, select variables, levels etc.
1236 if(writelevel(level)) then
1237 ! only output a grid when fully within clipped region selected
1238 ! by writespshift array
1239 if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1240 *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1241 <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
1242 morton_aim_p(morton_no)=.true.
1243 end if
1244 end if
1245 end do
1246 call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
1247 icomm,ierrmpi)
1248 select case(convert_type)
1249 case('vtuB64','vtuBmpi64')
1250 cell_corner=.true.
1251 case('vtuBCC64','vtuBCCmpi64')
1252 cell_corner=.false.
1253 end select
1254 if (mype /= 0) then
1255 do morton_no=morton_start(mype),morton_stop(mype)
1256 if(.not. morton_aim(morton_no)) cycle
1257 igrid=sfc_to_igrid(morton_no)
1258 call calc_x(igrid,xc,xcc)
1259 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1260 ixc^l,ixcc^l,.true.)
1261 itag=morton_no
1262 call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
1263 if(cell_corner) then
1264 call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
1265 else
1266 call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
1267 end if
1268 end do
1269 else
1270 ! mype==0
1271 offset=0
1272 inquire(qunit,opened=fileopen)
1273 if(.not.fileopen)then
1274 ! generate filename
1275 filenr=snapshotini
1276 if (autoconvert) filenr=snapshotnext
1277 write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
1278 ! Open the file for the header part
1279 open(qunit,file=filename,status='replace')
1280 end if
1281 call getheadernames(wnamei,xandwnamei,outfilehead)
1282 ! generate xml header
1283 write(qunit,'(a)')'<?xml version="1.0"?>'
1284 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
1285 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1286 write(qunit,'(a)')'<UnstructuredGrid>'
1287 write(qunit,'(a)')'<FieldData>'
1288 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1289 'NumberOfTuples="1" format="ascii">'
1290 write(qunit,*) real(global_time*time_convert_factor)
1291 write(qunit,'(a)')'</DataArray>'
1292 write(qunit,'(a)')'</FieldData>'
1293 ! number of cells, number of corner points, per grid.
1294 nx^d=ixmhi^d-ixmlo^d+1;
1295 nxc^d=nx^d+1;
1296 nc={nx^d*}
1297 np={nxc^d*}
1298 length=np*size_double
1299 lengthcc=nc*size_double
1300 length_coords=3*length
1301 length_conn=2**^nd*size_int*nc
1302 length_offsets=nc*size_int
1303 ! Note: using the w_write, writelevel, writespshift
1304 do morton_no=morton_start(0),morton_stop(0)
1305 if(.not. morton_aim(morton_no)) cycle
1306 if(cell_corner) then
1307 ! we write out every grid as one VTK PIECE
1308 write(qunit,'(a,i7,a,i7,a)') &
1309 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1310 write(qunit,'(a)')'<PointData>'
1311 do iw=1,nw+nwauxio
1312 if(iw<=nw) then
1313 if(.not.w_write(iw)) cycle
1314 end if
1315 write(qunit,'(a,a,a,i16,a)')&
1316 '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
1317 '" format="appended" offset="',offset,'">'
1318 write(qunit,'(a)')'</DataArray>'
1319 offset=offset+length+size_int
1320 end do
1321 write(qunit,'(a)')'</PointData>'
1322 write(qunit,'(a)')'<Points>'
1323 write(qunit,'(a,i16,a)') &
1324 '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1325 ! write cell corner coordinates in a backward dimensional loop, always 3D output
1326 offset=offset+length_coords+size_int
1327 write(qunit,'(a)')'</Points>'
1328 else
1329 ! we write out every grid as one VTK PIECE
1330 write(qunit,'(a,i7,a,i7,a)') &
1331 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1332 write(qunit,'(a)')'<CellData>'
1333 do iw=1,nw+nwauxio
1334 if(iw<=nw) then
1335 if(.not.w_write(iw)) cycle
1336 end if
1337 write(qunit,'(a,a,a,i16,a)')&
1338 '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
1339 '" format="appended" offset="',offset,'">'
1340 write(qunit,'(a)')'</DataArray>'
1341 offset=offset+lengthcc+size_int
1342 end do
1343 write(qunit,'(a)')'</CellData>'
1344 write(qunit,'(a)')'<Points>'
1345 write(qunit,'(a,i16,a)') &
1346 '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1347 ! write cell corner coordinates in a backward dimensional loop, always 3D output
1348 offset=offset+length_coords+size_int
1349 write(qunit,'(a)')'</Points>'
1350 end if
1351 write(qunit,'(a)')'<Cells>'
1352 ! connectivity part
1353 write(qunit,'(a,i16,a)')&
1354 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
1355 offset=offset+length_conn+size_int
1356 ! offsets data array
1357 write(qunit,'(a,i16,a)') &
1358 '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
1359 offset=offset+length_offsets+size_int
1360 ! VTK cell type data array
1361 write(qunit,'(a,i16,a)') &
1362 '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
1363 offset=offset+size_int+nc*size_int
1364 write(qunit,'(a)')'</Cells>'
1365 write(qunit,'(a)')'</Piece>'
1366 end do
1367 ! write metadata communicated from other processors
1368 if(npe>1)then
1369 do ipe=1, npe-1
1370 do morton_no=morton_start(ipe),morton_stop(ipe)
1371 if(.not. morton_aim(morton_no)) cycle
1372 if(cell_corner) then
1373 ! we write out every grid as one VTK PIECE
1374 write(qunit,'(a,i7,a,i7,a)') &
1375 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1376 write(qunit,'(a)')'<PointData>'
1377 do iw=1,nw+nwauxio
1378 if(iw<=nw) then
1379 if(.not.w_write(iw)) cycle
1380 end if
1381 write(qunit,'(a,a,a,i16,a)')&
1382 '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
1383 '" format="appended" offset="',offset,'">'
1384 write(qunit,'(a)')'</DataArray>'
1385 offset=offset+length+size_int
1386 end do
1387 write(qunit,'(a)')'</PointData>'
1388 write(qunit,'(a)')'<Points>'
1389 write(qunit,'(a,i16,a)') &
1390 '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1391 ! write cell corner coordinates in a backward dimensional loop, always 3D output
1392 offset=offset+length_coords+size_int
1393 write(qunit,'(a)')'</Points>'
1394 else
1395 ! we write out every grid as one VTK PIECE
1396 write(qunit,'(a,i7,a,i7,a)') &
1397 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1398 write(qunit,'(a)')'<CellData>'
1399 do iw=1,nw+nwauxio
1400 if(iw<=nw) then
1401 if(.not.w_write(iw)) cycle
1402 end if
1403 write(qunit,'(a,a,a,i16,a)')&
1404 '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
1405 '" format="appended" offset="',offset,'">'
1406 write(qunit,'(a)')'</DataArray>'
1407 offset=offset+lengthcc+size_int
1408 end do
1409 write(qunit,'(a)')'</CellData>'
1410 write(qunit,'(a)')'<Points>'
1411 write(qunit,'(a,i16,a)') &
1412 '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1413 ! write cell corner coordinates in a backward dimensional loop, always 3D output
1414 offset=offset+length_coords+size_int
1415 write(qunit,'(a)')'</Points>'
1416 end if
1417 write(qunit,'(a)')'<Cells>'
1418 ! connectivity part
1419 write(qunit,'(a,i16,a)')&
1420 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
1421 offset=offset+length_conn+size_int
1422 ! offsets data array
1423 write(qunit,'(a,i16,a)') &
1424 '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
1425 offset=offset+length_offsets+size_int
1426 ! VTK cell type data array
1427 write(qunit,'(a,i16,a)') &
1428 '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
1429 offset=offset+size_int+nc*size_int
1430 write(qunit,'(a)')'</Cells>'
1431 write(qunit,'(a)')'</Piece>'
1432 end do
1433 end do
1434 end if
1435 write(qunit,'(a)')'</UnstructuredGrid>'
1436 write(qunit,'(a)')'<AppendedData encoding="raw">'
1437 close(qunit)
1438 open(qunit,file=filename,access='stream',form='unformatted',position='append')
1439 buf='_'
1440 write(qunit) trim(buf)
1441 do morton_no=morton_start(0),morton_stop(0)
1442 if(.not. morton_aim(morton_no)) cycle
1443 igrid=sfc_to_igrid(morton_no)
1444 call calc_x(igrid,xc,xcc)
1445 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1446 ixc^l,ixcc^l,.true.)
1447 do iw=1,nw+nwauxio
1448 if(iw<=nw) then
1449 if(.not.w_write(iw)) cycle
1450 end if
1451 if(cell_corner) then
1452 write(qunit) length
1453 write(qunit) {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
1454 else
1455 write(qunit) lengthcc
1456 write(qunit) {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
1457 end if
1458 end do
1459 write(qunit) length_coords
1460 {do ix^db=ixcmin^db,ixcmax^db \}
1461 x_vtk(1:3)=zero;
1462 x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1463 do k=1,3
1464 write(qunit) x_vtk(k)
1465 end do
1466 {end do \}
1467 write(qunit) length_conn
1468 {do ix^db=1,nx^db\}
1469 {^ifoned write(qunit)ix1-1,ix1 \}
1470 {^iftwod
1471 write(qunit)(ix2-1)*nxc1+ix1-1, &
1472 (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1473 \}
1474 {^ifthreed
1475 write(qunit)&
1476 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1477 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1478 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1479 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1480 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1481 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1482 ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1483 ix3*nxc2*nxc1+ ix2*nxc1+ix1
1484 \}
1485 {end do\}
1486 write(qunit) length_offsets
1487 do icel=1,nc
1488 write(qunit) icel*(2**^nd)
1489 end do
1490 {^ifoned vtk_type=3 \}
1491 {^iftwod vtk_type=8 \}
1492 {^ifthreed vtk_type=11 \}
1493 write(qunit) size_int*nc
1494 do icel=1,nc
1495 write(qunit) vtk_type
1496 end do
1497 end do
1498 allocate(intstatus(mpi_status_size,1))
1499 if(npe>1)then
1500 ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d;
1501 ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
1502 do ipe=1, npe-1
1503 do morton_no=morton_start(ipe),morton_stop(ipe)
1504 if(.not. morton_aim(morton_no)) cycle
1505 itag=morton_no
1506 call mpi_recv(xc_tmp,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1507 if(cell_corner) then
1508 call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1509 else
1510 call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1511 end if
1512 do iw=1,nw+nwauxio
1513 if(iw<=nw) then
1514 if(.not.w_write(iw)) cycle
1515 end if
1516 if(cell_corner) then
1517 write(qunit) length
1518 write(qunit) {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
1519 else
1520 write(qunit) lengthcc
1521 write(qunit) {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
1522 end if
1523 end do
1524 write(qunit) length_coords
1525 {do ix^db=ixcmin^db,ixcmax^db \}
1526 x_vtk(1:3)=zero;
1527 x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1528 do k=1,3
1529 write(qunit) x_vtk(k)
1530 end do
1531 {end do \}
1532 write(qunit) length_conn
1533 {do ix^db=1,nx^db\}
1534 {^ifoned write(qunit)ix1-1,ix1 \}
1535 {^iftwod
1536 write(qunit)(ix2-1)*nxc1+ix1-1, &
1537 (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1538 \}
1539 {^ifthreed
1540 write(qunit)&
1541 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1542 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1543 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1544 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1545 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1546 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1547 ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1548 ix3*nxc2*nxc1+ ix2*nxc1+ix1
1549 \}
1550 {end do\}
1551 write(qunit) length_offsets
1552 do icel=1,nc
1553 write(qunit) icel*(2**^nd)
1554 end do
1555 {^ifoned vtk_type=3 \}
1556 {^iftwod vtk_type=8 \}
1557 {^ifthreed vtk_type=11 \}
1558 write(qunit) size_int*nc
1559 do icel=1,nc
1560 write(qunit) vtk_type
1561 end do
1562 end do
1563 end do
1564 end if
1565 close(qunit)
1566 open(qunit,file=filename,status='unknown',form='formatted',position='append')
1567 write(qunit,'(a)')'</AppendedData>'
1568 write(qunit,'(a)')'</VTKFile>'
1569 close(qunit)
1570 deallocate(intstatus)
1571 end if
1572 deallocate(morton_aim,morton_aim_p)
1573 if (npe>1) then
1574 call mpi_barrier(icomm,ierrmpi)
1575 end if
1576
1577 end subroutine unstructuredvtkb64
1578
1579 subroutine save_connvtk(qunit,igrid)
1580 ! this saves the basic line, pixel and voxel connectivity,
1581 ! as used by VTK file outputs for unstructured grid
1583
1584 integer, intent(in) :: qunit, igrid
1585
1586 integer :: nx^D, nxC^D, ix^D
1587
1588 nx^d=ixmhi^d-ixmlo^d+1;
1589 nxc^d=nx^d+1;
1590 {do ix^db=1,nx^db\}
1591 {^ifoned write(qunit,'(2(i7,1x))')ix1-1,ix1 \}
1592 {^iftwod
1593 write(qunit,'(4(i7,1x))')(ix2-1)*nxc1+ix1-1, &
1594 (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1595 \}
1596 {^ifthreed
1597 write(qunit,'(8(i7,1x))')&
1598 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1599 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1600 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1601 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1602 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1603 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1604 ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1605 ix3*nxc2*nxc1+ ix2*nxc1+ix1
1606 \}
1607 {end do\}
1608
1609 end subroutine save_connvtk
1610
1611 subroutine imagedatavtk_mpi(qunit)
1612 ! output for vti format to paraview, non-binary version output
1613 ! parallel, uses calc_grid to compute nwauxio variables
1614 ! allows renormalizing using convert factors
1615 ! allows skipping of w_write selected variables
1616 ! implementation such that length of ASCII output is identical when
1617 ! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
1621
1622 integer, intent(in) :: qunit
1623
1624 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
1625 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
1626 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1627 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1628 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
1629 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
1630 double precision, dimension(0:nw+nwauxio) :: normconv
1631 double precision :: origin(1:3), spacing(1:3)
1632 integer :: igrid,iigrid,level,ixC^L,ixCC^L
1633 integer :: NumGridsOnLevel(1:nlevelshi)
1634 integer :: nx^D
1635 integer :: filenr
1636 integer :: itag,ipe,Morton_no,Morton_length
1637 integer :: ixrvC^L, ixrvCC^L, siz_ind, ind_send(5*^ND), ind_recv(5*^ND)
1638 integer :: wholeExtent(1:6), ig^D
1639 integer, allocatable :: intstatus(:,:)
1640 logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
1641 logical :: fileopen
1642 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1643 character(len=1024) :: outfilehead
1644 character(len=80):: filename
1645 type(tree_node_ptr) :: tree
1646
1647 if(levmin/=levmax) call mpistop('ImageData can only be used when levmin=levmax')
1648 normconv(0) = length_convert_factor
1649 normconv(1:nw) = w_convert_factor
1650 siz_ind=5*^nd
1651 morton_length=morton_stop(npe-1)-morton_start(0)+1
1652 allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
1653 allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
1654 morton_aim=.false.
1655 morton_aim_p=.false.
1656 do morton_no=morton_start(mype),morton_stop(mype)
1657 igrid=sfc_to_igrid(morton_no)
1658 level=node(plevel_,igrid)
1659 ! we can clip parts of the grid away, select variables, levels etc.
1660 if(writelevel(level)) then
1661 ! only output a grid when fully within clipped region selected
1662 ! by writespshift array
1663 if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1664 *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1665 <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
1666 morton_aim_p(morton_no)=.true.
1667 end if
1668 end if
1669 end do
1670 call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
1671 icomm,ierrmpi)
1672 if(mype /= 0) then
1673 do morton_no=morton_start(mype),morton_stop(mype)
1674 if(.not. morton_aim(morton_no)) cycle
1675 igrid=sfc_to_igrid(morton_no)
1676 call calc_x(igrid,xc,xcc)
1677 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1678 ixc^l,ixcc^l,.true.)
1679 tree%node => igrid_to_node(igrid, mype)%node
1680 {^d& ig^d = tree%node%ig^d; }
1681 itag=morton_no
1682 ind_send=(/ ixc^l,ixcc^l, ig^d /)
1683 call mpi_send(ind_send,siz_ind,mpi_integer, 0,itag,icomm,ierrmpi)
1684 call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
1685 call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
1686 end do
1687 else
1688 inquire(qunit,opened=fileopen)
1689 if(.not.fileopen)then
1690 ! generate filename
1691 filenr=snapshotini
1692 if (autoconvert) filenr=snapshotnext
1693 write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vti"
1694 ! Open the file for the header part
1695 open(qunit,file=filename,status='unknown',form='formatted')
1696 end if
1697 call getheadernames(wnamei,xandwnamei,outfilehead)
1698 ! number of cells per grid.
1699 nx^d=ixmhi^d-ixmlo^d+1;
1700 origin = 0
1701 {^d& origin(^d) = xprobmin^d*normconv(0); }
1702 spacing = zero
1703 {^d&spacing(^d) = dxlevel(^d)*normconv(0); }
1704 wholeextent = 0
1705 ! if we use writespshift, the whole extent has to be calculated:
1706 {^d&wholeextent(^d*2-1) = nx^d * ceiling(((xprobmax^d-xprobmin^d)*writespshift(^d,1)) &
1707 /(nx^d*dxlevel(^d))) \}
1708 {^d&wholeextent(^d*2) = nx^d * floor(((xprobmax^d-xprobmin^d)*(1.0d0-writespshift(^d,2))) &
1709 /(nx^d*dxlevel(^d))) \}
1710
1711 ! generate xml header
1712 write(qunit,'(a)')'<?xml version="1.0"?>'
1713 write(qunit,'(a)',advance='no') '<VTKFile type="ImageData"'
1714 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1715 write(qunit,'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')' <ImageData Origin="',&
1716 origin,'" WholeExtent="',wholeextent,'" Spacing="',spacing,'">'
1717 write(qunit,'(a)')'<FieldData>'
1718 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1719 'NumberOfTuples="1" format="ascii">'
1720 write(qunit,*) real(global_time*time_convert_factor)
1721 write(qunit,'(a)')'</DataArray>'
1722 write(qunit,'(a)')'</FieldData>'
1723
1724 ! write the data from proc 0
1725 do morton_no=morton_start(0),morton_stop(0)
1726 if(.not. morton_aim(morton_no)) cycle
1727 igrid=sfc_to_igrid(morton_no)
1728 tree%node => igrid_to_node(igrid, 0)%node
1729 {^d& ig^d = tree%node%ig^d; }
1730 call calc_x(igrid,xc,xcc)
1731 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1732 ixc^l,ixcc^l,.true.)
1733 call write_vti(qunit,ixg^ll,ixc^l,ixcc^l,ig^d,&
1734 nx^d,normconv,wnamei,wc_tmp,wcc_tmp)
1735 end do
1736
1737 if(npe>1)then
1738 allocate(intstatus(mpi_status_size,1))
1739 do ipe=1, npe-1
1740 do morton_no=morton_start(ipe),morton_stop(ipe)
1741 if(.not. morton_aim(morton_no)) cycle
1742 itag=morton_no
1743 call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1744 ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
1745 ixrvccmin^d=ind_recv(2*^nd+^d);ixrvccmax^d=ind_recv(3*^nd+^d);
1746 ig^d=ind_recv(4*^nd+^d);
1747 call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1748 call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1749 call write_vti(qunit,ixg^ll,ixrvc^l,ixrvcc^l,ig^d,&
1750 nx^d,normconv,wnamei,wc_tmp,wcc_tmp)
1751 end do
1752 end do
1753 end if
1754 write(qunit,'(a)')'</ImageData>'
1755 write(qunit,'(a)')'</VTKFile>'
1756 close(qunit)
1757 if(npe>1) deallocate(intstatus)
1758 end if
1759
1760 deallocate(morton_aim,morton_aim_p)
1761 if (npe>1) then
1762 call mpi_barrier(icomm,ierrmpi)
1763 endif
1764
1765 end subroutine imagedatavtk_mpi
1766
1767 subroutine punstructuredvtk_mpi(qunit)
1768 ! Write one pvtu and vtu files for each processor
1769 ! Otherwise like unstructuredvtk_mpi
1773
1774 integer, intent(in) :: qunit
1775
1776 double precision, dimension(0:nw+nwauxio) :: normconv
1777 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
1778 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
1779 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1780 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1781 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
1782 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
1783 integer :: nx^D,nxC^D,nc,np, igrid,ixC^L,ixCC^L,level,Morton_no
1784 integer :: filenr
1785 logical :: fileopen,conv_grid
1786 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1787 character(len=1024) :: outfilehead
1788 character(len=80) :: pfilename
1789
1790 ! Write pvtu-file:
1791 if (mype==0) then
1792 call write_pvtu(qunit)
1793 endif
1794 ! Now write the Source files:
1795 inquire(qunit,opened=fileopen)
1796 if(.not.fileopen)then
1797 ! generate filename
1798 filenr=snapshotini
1799 if (autoconvert) filenr=snapshotnext
1800 ! Open the file for the header part
1801 write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename),filenr,"p",mype,".vtu"
1802 open(qunit,file=pfilename,status='unknown',form='formatted')
1803 end if
1804 ! generate xml header
1805 write(qunit,'(a)')'<?xml version="1.0"?>'
1806 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
1807 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1808 write(qunit,'(a)')' <UnstructuredGrid>'
1809 write(qunit,'(a)')'<FieldData>'
1810 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1811 'NumberOfTuples="1" format="ascii">'
1812 write(qunit,*) real(global_time*time_convert_factor)
1813 write(qunit,'(a)')'</DataArray>'
1814 write(qunit,'(a)')'</FieldData>'
1815
1816 call getheadernames(wnamei,xandwnamei,outfilehead)
1817
1818 ! number of cells, number of corner points, per grid.
1819 nx^d=ixmhi^d-ixmlo^d+1;
1820 nxc^d=nx^d+1;
1821 nc={nx^d*}
1822 np={nxc^d*}
1823
1824 ! Note: using the w_write, writelevel, writespshift
1825 ! we can clip parts of the grid away, select variables, levels etc.
1826 do level=levmin,levmax
1827 if (.not.writelevel(level)) cycle
1828 do morton_no=morton_start(mype),morton_stop(mype)
1829 igrid=sfc_to_igrid(morton_no)
1830 if (node(plevel_,igrid)/=level) cycle
1831 ! only output a grid when fully within clipped region selected
1832 ! by writespshift array
1833 conv_grid=({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1834 *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1835 <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})
1836 if (.not.conv_grid) cycle
1837 call calc_x(igrid,xc,xcc)
1838 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1839 ixc^l,ixcc^l,.true.)
1840 call write_vtk(qunit,ixg^ll,ixc^l,ixcc^l,igrid,nc,np,nx^d,nxc^d,&
1841 normconv,wnamei,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp)
1842 end do ! Morton_no loop
1843 end do ! level loop
1844
1845 write(qunit,'(a)')' </UnstructuredGrid>'
1846 write(qunit,'(a)')'</VTKFile>'
1847 close(qunit)
1848
1849 if (npe>1) then
1850 call mpi_barrier(icomm,ierrmpi)
1851 end if
1852
1853 end subroutine punstructuredvtk_mpi
1854
1855 subroutine unstructuredvtk_mpi(qunit)
1856 ! output for vtu format to paraview, non-binary version output
1857 ! parallel, uses calc_grid to compute nwauxio variables
1858 ! allows renormalizing using convert factors
1859 ! allows skipping of w_write selected variables
1860 ! implementation such that length of ASCII output is identical when
1861 ! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
1865
1866 integer, intent(in) :: qunit
1867
1868 double precision :: x_VTK(1:3)
1869 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
1870 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
1871 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1872 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1873 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
1874 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
1875 double precision, dimension(0:nw+nwauxio) :: normconv
1876 integer:: igrid,iigrid,level,ixC^L,ixCC^L
1877 integer:: NumGridsOnLevel(1:nlevelshi)
1878 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,ix^D
1879 integer :: filenr
1880 integer :: itag,ipe,Morton_no,siz_ind
1881 integer :: ind_send(4*^ND),ind_recv(4*^ND)
1882 integer :: levmin_recv,levmax_recv,level_recv,igrid_recv,ixrvC^L,ixrvCC^L
1883 integer, allocatable :: intstatus(:,:)
1884 logical :: fileopen,conv_grid,cond_grid_recv
1885 character(len=80):: filename
1886 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1887 character(len=1024) :: outfilehead
1888
1889 if(mype==0) then
1890 inquire(qunit,opened=fileopen)
1891 if(.not.fileopen)then
1892 ! generate filename
1893 filenr=snapshotini
1894 if (autoconvert) filenr=snapshotnext
1895 write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
1896 ! Open the file for the header part
1897 open(qunit,file=filename,status='unknown',form='formatted')
1898 end if
1899 ! generate xml header
1900 write(qunit,'(a)')'<?xml version="1.0"?>'
1901 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
1902 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1903 write(qunit,'(a)')'<UnstructuredGrid>'
1904 write(qunit,'(a)')'<FieldData>'
1905 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1906 'NumberOfTuples="1" format="ascii">'
1907 write(qunit,*) real(global_time*time_convert_factor)
1908 write(qunit,'(a)')'</DataArray>'
1909 write(qunit,'(a)')'</FieldData>'
1910 end if
1911
1912 call getheadernames(wnamei,xandwnamei,outfilehead)
1913 ! number of cells, number of corner points, per grid.
1914 nx^d=ixmhi^d-ixmlo^d+1;
1915 nxc^d=nx^d+1;
1916 nc={nx^d*}
1917 np={nxc^d*}
1918 ! all slave processors send their minmal/maximal levels
1919 if (mype/=0) then
1920 if (morton_stop(mype)==0) call mpistop("nultag")
1921 itag=1000*morton_stop(mype)
1922 !print *,'ype,itag for levmin=',mype,itag,levmin
1923 call mpi_send(levmin,1,mpi_integer, 0,itag,icomm,ierrmpi)
1924 itag=2000*morton_stop(mype)
1925 !print *,'mype,itag for levmax=',mype,itag,levmax
1926 call mpi_send(levmax,1,mpi_integer, 0,itag,icomm,ierrmpi)
1927 end if
1928 ! Note: using the w_write, writelevel, writespshift
1929 ! we can clip parts of the grid away, select variables, levels etc.
1930 do level=levmin,levmax
1931 if (.not.writelevel(level)) cycle
1932 do morton_no=morton_start(mype),morton_stop(mype)
1933 igrid=sfc_to_igrid(morton_no)
1934 if (mype/=0)then
1935 itag=morton_no
1936 call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
1937 itag=igrid
1938 call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
1939 end if
1940 if (node(plevel_,igrid)/=level) cycle
1941 ! only output a grid when fully within clipped region selected
1942 ! by writespshift array
1943 conv_grid=({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1944 *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1945 <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})
1946 if (mype/=0)then
1947 call mpi_send(conv_grid,1,mpi_logical,0,itag,icomm,ierrmpi)
1948 end if
1949 if (.not.conv_grid) cycle
1950 call calc_x(igrid,xc,xcc)
1951 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1952 ixc^l,ixcc^l,.true.)
1953 if(mype/=0) then
1954 itag=morton_no
1955 ind_send=(/ ixc^l,ixcc^l /)
1956 siz_ind=4*^nd
1957 call mpi_send(ind_send,siz_ind,mpi_integer, 0,itag,icomm,ierrmpi)
1958 call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
1959 call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
1960 call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
1961 itag=igrid
1962 call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
1963 call mpi_send(xcc_tmp,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
1964 else
1965 call write_vtk(qunit,ixg^ll,ixc^l,ixcc^l,igrid,nc,np,nx^d,nxc^d,&
1966 normconv,wnamei,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp)
1967 end if
1968 end do ! Morton_no loop
1969 end do ! level loop
1970
1971 if(mype==0) then
1972 allocate(intstatus(mpi_status_size,1))
1973 if(npe>1)then
1974 do ipe=1,npe-1
1975 itag=1000*morton_stop(ipe)
1976 call mpi_recv(levmin_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1977 !!print *,'mype RECEIVES,itag for levmin=',mype,itag,levmin_recv
1978 itag=2000*morton_stop(ipe)
1979 call mpi_recv(levmax_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1980 !!print *,'mype RECEIVES itag for levmax=',mype,itag,levmax_recv
1981 do level=levmin_recv,levmax_recv
1982 if (.not.writelevel(level)) cycle
1983 do morton_no=morton_start(ipe),morton_stop(ipe)
1984 itag=morton_no
1985 call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1986 itag=igrid_recv
1987 call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1988 if (level_recv/=level) cycle
1989 call mpi_recv(cond_grid_recv,1,mpi_logical, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1990 if(.not.cond_grid_recv)cycle
1991 itag=morton_no
1992 siz_ind=4*^nd
1993 call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1994 ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
1995 ixrvccmin^d=ind_recv(2*^nd+^d);ixrvccmax^d=ind_recv(3*^nd+^d);
1996 call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
1997 ,icomm,intstatus(:,1),ierrmpi)
1998 call mpi_recv(wc_tmp_recv,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1999 call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2000 itag=igrid_recv
2001 call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2002 call mpi_recv(xcc_tmp_recv,1,type_block_xcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2003 call write_vtk(qunit,ixg^ll,ixrvc^l,ixrvcc^l,igrid_recv,&
2004 nc,np,nx^d,nxc^d,normconv,wnamei,&
2005 xc_tmp_recv,xcc_tmp_recv,wc_tmp_recv,wcc_tmp_recv)
2006 end do ! Morton_no loop
2007 end do ! level loop
2008 end do ! processor loop
2009 end if ! multiple processors
2010 write(qunit,'(a)')'</UnstructuredGrid>'
2011 write(qunit,'(a)')'</VTKFile>'
2012 close(qunit)
2013 end if
2014 if (npe>1) then
2015 call mpi_barrier(icomm,ierrmpi)
2016 if(mype==0)deallocate(intstatus)
2017 end if
2018
2019 end subroutine unstructuredvtk_mpi
2020
2021 subroutine write_vtk(qunit,ixI^L,ixC^L,ixCC^L,igrid,nc,np,nx^D,nxC^D,&
2022 normconv,wnamei,xC,xCC,wC,wCC)
2024
2025 integer, intent(in) :: qunit
2026 integer, intent(in) :: ixI^L,ixC^L,ixCC^L
2027 integer, intent(in) :: igrid,nc,np,nx^D,nxC^D
2028 double precision, intent(in) :: normconv(0:nw+nwauxio)
2029 character(len=name_len), intent(in):: wnamei(1:nw+nwauxio)
2030 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
2031 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
2032 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC
2033 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC
2034
2035 double precision :: x_VTK(1:3)
2036 integer :: iw,ix^D,icel,VTK_type
2037
2038 select case(convert_type)
2039 case('vtumpi','pvtumpi')
2040 ! we write out every grid as one VTK PIECE
2041 write(qunit,'(a,i7,a,i7,a)') &
2042 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2043 write(qunit,'(a)')'<PointData>'
2044 do iw=1,nw+nwauxio
2045 if(iw<=nw) then
2046 if(.not.w_write(iw)) cycle
2047 end if
2048 write(qunit,'(a,a,a)')&
2049 '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
2050 write(qunit,'(200(1pe14.6))') {(|}wc(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
2051 write(qunit,'(a)')'</DataArray>'
2052 end do
2053 write(qunit,'(a)')'</PointData>'
2054 write(qunit,'(a)')'<Points>'
2055 write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
2056 ! write cell corner coordinates in a backward dimensional loop, always 3D output
2057 {do ix^db=ixcmin^db,ixcmax^db \}
2058 x_vtk(1:3)=zero;
2059 x_vtk(1:ndim)=xc(ix^d,1:ndim)*normconv(0);
2060 write(qunit,'(3(1pe14.6))') x_vtk
2061 {end do \}
2062 write(qunit,'(a)')'</DataArray>'
2063 write(qunit,'(a)')'</Points>'
2064
2065 case('vtuCCmpi','pvtuCCmpi')
2066 ! we write out every grid as one VTK PIECE
2067 write(qunit,'(a,i7,a,i7,a)') &
2068 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2069 write(qunit,'(a)')'<CellData>'
2070 do iw=1,nw+nwauxio
2071 if(iw<=nw) then
2072 if(.not.w_write(iw)) cycle
2073 end if
2074 write(qunit,'(a,a,a)')&
2075 '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
2076 write(qunit,'(200(1pe14.6))') {(|}wcc(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
2077 write(qunit,'(a)')'</DataArray>'
2078 end do
2079 write(qunit,'(a)')'</CellData>'
2080 write(qunit,'(a)')'<Points>'
2081 write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
2082 ! write cell corner coordinates in a backward dimensional loop, always 3D output
2083 {do ix^db=ixcmin^db,ixcmax^db \}
2084 x_vtk(1:3)=zero;
2085 x_vtk(1:ndim)=xc(ix^d,1:ndim)*normconv(0);
2086 write(qunit,'(3(1pe14.6))') x_vtk
2087 {end do \}
2088 write(qunit,'(a)')'</DataArray>'
2089 write(qunit,'(a)')'</Points>'
2090 end select
2091
2092 write(qunit,'(a)')'<Cells>'
2093 ! connectivity part
2094 write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
2095 call save_connvtk(qunit,igrid)
2096 write(qunit,'(a)')'</DataArray>'
2097 ! offsets data array
2098 write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
2099 do icel=1,nc
2100 write(qunit,'(i7)') icel*(2**^nd)
2101 end do
2102 write(qunit,'(a)')'</DataArray>'
2103 ! VTK cell type data array
2104 write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
2105 ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
2106 {^ifoned vtk_type=3 \}
2107 {^iftwod vtk_type=8 \}
2108 {^ifthreed vtk_type=11 \}
2109 do icel=1,nc
2110 write(qunit,'(i2)') vtk_type
2111 end do
2112 write(qunit,'(a)')'</DataArray>'
2113 write(qunit,'(a)')'</Cells>'
2114 write(qunit,'(a)')'</Piece>'
2115
2116 end subroutine write_vtk
2117
2118 subroutine write_vti(qunit,ixI^L,ixC^L,ixCC^L,ig^D,nx^D,&
2119 normconv,wnamei,wC,wCC)
2121
2122 integer, intent(in) :: qunit
2123 integer, intent(in) :: ixI^L,ixC^L,ixCC^L
2124 integer, intent(in) :: ig^D,nx^D
2125 double precision, intent(in) :: normconv(0:nw+nwauxio)
2126 character(len=name_len), intent(in):: wnamei(1:nw+nwauxio)
2127 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC
2128 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC
2129
2130 integer :: iw,ix^D
2131 integer :: extent(1:6)
2132
2133 extent = 0
2134 {^d& extent(^d*2-1) = (ig^d-1) * nx^d; }
2135 {^d& extent(^d*2) = (ig^d) * nx^d; }
2136
2137 select case(convert_type)
2138 case('vtimpi','pvtimpi')
2139 ! we write out every grid as one VTK PIECE
2140 write(qunit,'(a,6(i10),a)') &
2141 '<Piece Extent="',extent,'">'
2142 write(qunit,'(a)')'<PointData>'
2143 do iw=1,nw+nwauxio
2144 if(iw<=nw) then
2145 if(.not.w_write(iw)) cycle
2146 end if
2147 write(qunit,'(a,a,a)')&
2148 '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
2149 write(qunit,'(200(1pe20.12))') {(|}wc(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
2150 write(qunit,'(a)')'</DataArray>'
2151 end do
2152 write(qunit,'(a)')'</PointData>'
2153 case('vtiCCmpi','pvtiCCmpi')
2154 ! we write out every grid as one VTK PIECE
2155 write(qunit,'(a,6(i10),a)') &
2156 '<Piece Extent="',extent,'">'
2157 write(qunit,'(a)')'<CellData>'
2158 do iw=1,nw+nwauxio
2159 if(iw<=nw) then
2160 if(.not.w_write(iw)) cycle
2161 end if
2162 write(qunit,'(a,a,a)')&
2163 '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
2164 write(qunit,'(200(1pe20.12))') {(|}wcc(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
2165 write(qunit,'(a)')'</DataArray>'
2166 end do
2167 write(qunit,'(a)')'</CellData>'
2168 end select
2169
2170 write(qunit,'(a)')'</Piece>'
2171
2172 end subroutine write_vti
2173
2174 subroutine write_pvtu(qunit)
2177
2178 integer, intent(in) :: qunit
2179
2180 integer :: filenr,iw,ipe,iscalars
2181 logical :: fileopen
2182 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio),outtype
2183 character(len=1024) :: outfilehead
2184 character(len=80) :: filename,pfilename
2185
2186 select case(convert_type)
2187 case('pvtumpi','pvtuBmpi')
2188 outtype="PPointData"
2189 case('pvtuCCmpi','pvtuBCCmpi')
2190 outtype="PCellData"
2191 end select
2192 inquire(qunit,opened=fileopen)
2193 if(.not.fileopen)then
2194 ! generate filename
2195 filenr=snapshotini
2196 if (autoconvert) filenr=snapshotnext
2197 write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".pvtu"
2198 ! Open the file
2199 open(qunit,file=filename,status='unknown',form='formatted')
2200 end if
2201
2202 call getheadernames(wnamei,xandwnamei,outfilehead)
2203 ! Get the default selection:
2204 iscalars=1
2205 do iw=nw,1, -1
2206 if (w_write(iw)) iscalars=iw
2207 end do
2208 ! generate xml header
2209 write(qunit,'(a)')'<?xml version="1.0"?>'
2210 write(qunit,'(a)',advance='no') '<VTKFile type="PUnstructuredGrid"'
2211 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
2212 write(qunit,'(a)')' <PUnstructuredGrid GhostLevel="0">'
2213 ! Either celldata or pointdata:
2214 write(qunit,'(a,a,a,a,a)')&
2215 ' <',trim(outtype),' Scalars="',trim(wnamei(iscalars))//'">'
2216 do iw=1,nw
2217 if(.not.w_write(iw))cycle
2218 write(qunit,'(a,a,a)')&
2219 ' <PDataArray type="Float32" Name="',trim(wnamei(iw)),'"/>'
2220 end do
2221 do iw=nw+1,nw+nwauxio
2222 write(qunit,'(a,a,a)')&
2223 ' <PDataArray type="Float32" Name="',trim(wnamei(iw)),'"/>'
2224 end do
2225 write(qunit,'(a,a,a)')' </',trim(outtype),'>'
2226 write(qunit,'(a)')' <PPoints>'
2227 write(qunit,'(a)')' <PDataArray type="Float32" NumberOfComponents="3"/>'
2228 write(qunit,'(a)')' </PPoints>'
2229
2230 do ipe=0,npe-1
2231 write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename(&
2232 index(base_filename, '/', back = .true.)+1:&
2233 len(base_filename))),filenr,"p",&
2234 ipe,".vtu"
2235 write(qunit,'(a,a,a)')' <Piece Source="',trim(pfilename),'"/>'
2236 end do
2237 write(qunit,'(a)')' </PUnstructuredGrid>'
2238 write(qunit,'(a)')'</VTKFile>'
2239 close(qunit)
2240
2241 end subroutine write_pvtu
2242
2243 subroutine tecplot_mpi(qunit)
2244 ! output for tecplot (ASCII format)
2245 ! parallel, uses calc_grid to compute nwauxio variables
2246 ! allows renormalizing using convert factors
2247 ! the current implementation is such that tecplotmpi and tecplotCCmpi will
2248 ! create different length output ASCII files when used on 1 versus multiple CPUs
2249 ! in fact, on 1 CPU, there will be as many zones as there are levels
2250 ! on multiple CPUs, there will be a number of zones up to the number of
2251 ! levels times the number of CPUs (can be less, when some level not on a CPU)
2255
2256 integer, intent(in) :: qunit
2257
2258 double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
2259 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
2260 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
2261 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
2262 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
2263 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
2264 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
2265 double precision, dimension(0:nw+nwauxio) :: normconv
2266 integer:: igrid,iigrid,level,igonlevel,iw,idim,ix^D
2267 integer:: NumGridsOnLevel(1:nlevelshi)
2268 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,ixC^L,ixCC^L
2269 integer :: nodesonlevelmype,elemsonlevelmype
2270 integer :: nodes, elems
2271 integer, allocatable :: intstatus(:,:)
2272 integer :: itag,Morton_no,ipe,levmin_recv,levmax_recv,igrid_recv,level_recv
2273 integer :: ixrvC^L,ixrvCC^L
2274 integer :: ind_send(2*^ND),ind_recv(2*^ND),siz_ind,igonlevel_recv
2275 integer :: NumGridsOnLevel_mype(1:nlevelshi,0:npe-1)
2276 integer :: filenr
2277 logical :: fileopen,first
2278 character(len=80) :: filename
2279 character(len=1024) :: tecplothead
2280 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
2281 character(len=1024) :: outfilehead
2282
2283 if(nw/=count(w_write(1:nw)))then
2284 if(mype==0) print *,'tecplot_mpi does not use w_write=F'
2285 call mpistop('w_write, tecplot')
2286 end if
2287
2288 if(nocartesian)then
2289 if(mype==0) print *,'tecplot_mpi with nocartesian'
2290 end if
2291
2292 master_cpu_open : if (mype == 0) then
2293 inquire(qunit,opened=fileopen)
2294 if (.not.fileopen) then
2295 ! generate filename
2296 filenr=snapshotini
2297 if (autoconvert) filenr=snapshotnext
2298 write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".plt"
2299 open(qunit,file=filename,status='unknown')
2300 end if
2301 call getheadernames(wnamei,xandwnamei,outfilehead)
2302 write(tecplothead,'(a)') "VARIABLES = "//trim(outfilehead)
2303 write(qunit,'(a)') tecplothead(1:len_trim(tecplothead))
2304 end if master_cpu_open
2305
2306 ! determine overall number of grids per level, and the same info per CPU
2307 numgridsonlevel(1:nlevelshi)=0
2308 do level=levmin,levmax
2309 numgridsonlevel(level)=0
2310 do morton_no=morton_start(mype),morton_stop(mype)
2311 igrid = sfc_to_igrid(morton_no)
2312 if (node(plevel_,igrid)/=level) cycle
2313 numgridsonlevel(level)=numgridsonlevel(level)+1
2314 end do
2315 numgridsonlevel_mype(level,0:npe-1)=0
2316 numgridsonlevel_mype(level,mype) = numgridsonlevel(level)
2317 call mpi_allreduce(mpi_in_place,numgridsonlevel_mype(level,0:npe-1),npe,mpi_integer,&
2318 mpi_max,icomm,ierrmpi)
2319 call mpi_allreduce(mpi_in_place,numgridsonlevel(level),1,mpi_integer,mpi_sum, &
2320 icomm,ierrmpi)
2321 end do
2322
2323 nx^d=ixmhi^d-ixmlo^d+1;
2324 nxc^d=nx^d+1;
2325
2326 if(mype==0.and.npe>1) allocate(intstatus(mpi_status_size,1))
2327
2328 {^ifoned
2329 if(convert_type=='teclinempi') then
2330 nodes=0
2331 elems=0
2332 do level=levmin,levmax
2333 nodes=nodes + numgridsonlevel(level)*{nxc^d*}
2334 elems=elems + numgridsonlevel(level)*{nx^d*}
2335 end do
2336
2337 if (mype==0) write(qunit,"(a,i7,a,1pe12.5,a)") &
2338 'ZONE T="all levels", I=',elems, &
2339 ', SOLUTIONTIME=',global_time*time_convert_factor,', F=POINT'
2340
2341 igonlevel=0
2342 do morton_no=morton_start(mype),morton_stop(mype)
2343 igrid = sfc_to_igrid(morton_no)
2344 call calc_x(igrid,xc,xcc)
2345 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,ixc^l,ixcc^l,.true.)
2346 if (mype==0) then
2347 {do ix^db=ixccmin^db,ixccmax^db\}
2348 x_tec(1:ndim)=xcc_tmp(ix^d,1:ndim)*normconv(0)
2349 w_tec(1:nw+nwauxio)=wcc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2350 write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2351 {end do\}
2352 else if (mype/=0) then
2353 itag=morton_no
2354 call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2355 call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision,0,itag,icomm,ierrmpi)
2356 call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
2357 call mpi_send(xcc_tmp,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
2358 end if
2359 end do
2360 if(mype==0) then
2361 do ipe=1,npe-1
2362 do morton_no=morton_start(ipe),morton_stop(ipe)
2363 itag=morton_no
2364 call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2365 call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,&
2366 itag,icomm,intstatus(:,1),ierrmpi)
2367 call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,&
2368 icomm,intstatus(:,1),ierrmpi)
2369 call mpi_recv(xcc_tmp_recv,1,type_block_xcc_io, ipe,itag,&
2370 icomm,intstatus(:,1),ierrmpi)
2371 {do ix^db=ixccmin^db,ixccmax^db\}
2372 x_tec(1:ndim)=xcc_tmp_recv(ix^d,1:ndim)*normconv(0)
2373 w_tec(1:nw+nwauxio)=wcc_tmp_recv(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2374 write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2375 {end do\}
2376 end do
2377 end do
2378 close(qunit)
2379 end if
2380 else
2381 }
2382 if(mype/=0) then
2383 itag=1000*morton_stop(mype)
2384 call mpi_send(levmin,1,mpi_integer, 0,itag,icomm,ierrmpi)
2385 itag=2000*morton_stop(mype)
2386 call mpi_send(levmax,1,mpi_integer, 0,itag,icomm,ierrmpi)
2387 end if
2388
2389 do level=levmin,levmax
2390 nodesonlevelmype=numgridsonlevel_mype(level,mype)*{nxc^d*}
2391 elemsonlevelmype=numgridsonlevel_mype(level,mype)*{nx^d*}
2392 nodesonlevel=numgridsonlevel(level)*{nxc^d*}
2393 elemsonlevel=numgridsonlevel(level)*{nx^d*}
2394 ! for all tecplot variants coded up here, we let the TECPLOT ZONES coincide
2395 ! with the AMR grid LEVEL. Other options would be
2396 ! let each grid define a zone: inefficient for TECPLOT internal workings
2397 ! hence not implemented
2398 ! let entire octree define 1 zone: no difference in interpolation
2399 ! properties across TECPLOT zones detected as yet, hence not done
2400 select case(convert_type)
2401 case('tecplotmpi')
2402 ! in this option, we store the corner coordinates, as well as the corner
2403 ! values of all variables (obtained by averaging). This allows POINT packaging,
2404 ! and thus we can save full grid info by using one call to calc_grid
2405 if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2406 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
2407 'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2408 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
2409 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2410 do morton_no=morton_start(mype),morton_stop(mype)
2411 igrid = sfc_to_igrid(morton_no)
2412 if (mype/=0)then
2413 itag=morton_no
2414 call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2415 itag=igrid
2416 call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2417 end if
2418 if (node(plevel_,igrid)/=level) cycle
2419 call calc_x(igrid,xc,xcc)
2420 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2421 ixc^l,ixcc^l,.true.)
2422 if (mype/=0) then
2423 itag=morton_no
2424 ind_send=(/ ixc^l /)
2425 siz_ind=2*^nd
2426 call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2427 call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2428
2429 call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
2430 call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
2431 else
2432 {do ix^db=ixcmin^db,ixcmax^db\}
2433 x_tec(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0)
2434 w_tec(1:nw+nwauxio)=wc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2435 write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2436 {end do\}
2437 end if
2438 end do
2439 case('tecplotCCmpi')
2440 ! in this option, we store the corner coordinates, and the cell center
2441 ! values of all variables. Due to this mix of corner/cell center, we must
2442 ! use BLOCK packaging, and thus we have enormous overhead by using
2443 ! calc_grid repeatedly to merely fill values of cell corner coordinates
2444 ! and cell center values per dimension, per variable
2445 if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
2446 if(nw+nwauxio==1)then
2447 ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
2448 ! and just set [ndim+1]
2449 if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2450 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
2451 'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2452 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2453 ndim+1,']=CELLCENTERED), ZONETYPE=', &
2454 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2455 else
2456 if(ndim+nw+nwauxio<10) then
2457 ! difference only in length of integer format specification for ndim+nw+nwauxio
2458 if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2459 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
2460 'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2461 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2462 ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2463 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2464 else
2465 if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2466 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
2467 'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2468 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2469 ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2470 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2471 end if
2472 end if
2473
2474 do idim=1,ndim
2475 first=(idim==1)
2476 do morton_no=morton_start(mype),morton_stop(mype)
2477 igrid = sfc_to_igrid(morton_no)
2478 if (mype/=0)then
2479 itag=morton_no*idim
2480 call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2481 itag=igrid*idim
2482 call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2483 end if
2484 if (node(plevel_,igrid)/=level) cycle
2485 call calc_x(igrid,xc,xcc)
2486 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2487 ixc^l,ixcc^l,first)
2488 if (mype/=0)then
2489 ind_send=(/ ixc^l /)
2490 siz_ind=2*^nd
2491 itag=igrid*idim
2492 call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2493 call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2494 call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
2495 else
2496 write(qunit,fmt="(100(e14.6))") xc_tmp(ixc^s,idim)*normconv(0)
2497 end if
2498 end do
2499 end do
2500 do iw=1,nw+nwauxio
2501 do morton_no=morton_start(mype),morton_stop(mype)
2502 igrid = sfc_to_igrid(morton_no)
2503 if(mype/=0)then
2504 itag=morton_no*(ndim+iw)
2505 call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2506 itag=igrid*(ndim+iw)
2507 call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2508 end if
2509 if (node(plevel_,igrid)/=level) cycle
2510 call calc_x(igrid,xc,xcc)
2511 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2512 ixc^l,ixcc^l,.true.)
2513 if(mype/=0)then
2514 ind_send=(/ ixcc^l /)
2515 siz_ind=2*^nd
2516 itag=igrid*(ndim+iw)
2517 call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2518 call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2519 call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
2520 else
2521 write(qunit,fmt="(100(e14.6))") wcc_tmp(ixcc^s,iw)*normconv(iw)
2522 end if
2523 end do
2524 end do
2525 case default
2526 call mpistop('no such tecplot type')
2527 end select
2528
2529 igonlevel=0
2530 do morton_no=morton_start(mype),morton_stop(mype)
2531 igrid = sfc_to_igrid(morton_no)
2532 if(mype/=0)then
2533 itag=morton_no
2534 call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2535 itag=igrid
2536 call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2537 end if
2538 if(node(plevel_,igrid)/=level) cycle
2539 igonlevel=igonlevel+1
2540 if(mype/=0)then
2541 itag=igrid
2542 call mpi_send(igonlevel,1,mpi_integer, 0,itag,icomm,ierrmpi)
2543 end if
2544 if(mype==0)then
2545 call save_conntec(qunit,igrid,igonlevel)
2546 end if
2547 end do
2548 end do
2549
2550 if(mype==0 .and.npe>1) then
2551 do ipe=1,npe-1
2552 itag=1000*morton_stop(ipe)
2553 call mpi_recv(levmin_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2554 itag=2000*morton_stop(ipe)
2555 call mpi_recv(levmax_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2556 do level=levmin_recv,levmax_recv
2557 nodesonlevelmype=numgridsonlevel_mype(level,ipe)*{nxc^d*}
2558 elemsonlevelmype=numgridsonlevel_mype(level,ipe)*{nx^d*}
2559 nodesonlevel=numgridsonlevel(level)*{nxc^d*}
2560 elemsonlevel=numgridsonlevel(level)*{nx^d*}
2561 select case(convert_type)
2562 case('tecplotmpi')
2563 ! in this option, we store the corner coordinates, as well as the corner
2564 ! values of all variables (obtained by averaging). This allows POINT packaging,
2565 ! and thus we can save full grid info by using one call to calc_grid
2566 if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2567 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
2568 'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2569 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
2570 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2571 do morton_no=morton_start(ipe),morton_stop(ipe)
2572 itag=morton_no
2573 call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2574 itag=igrid_recv
2575 call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2576 if (level_recv/=level) cycle
2577 itag=morton_no
2578 siz_ind=2*^nd
2579 call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,&
2580 icomm,intstatus(:,1),ierrmpi)
2581 ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
2582 call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2583 ,icomm,intstatus(:,1),ierrmpi)
2584 call mpi_recv(wc_tmp_recv,1,type_block_wc_io, ipe,itag,&
2585 icomm,intstatus(:,1),ierrmpi)
2586 call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,&
2587 icomm,intstatus(:,1),ierrmpi)
2588 {do ix^db=ixrvcmin^db,ixrvcmax^db\}
2589 x_tec(1:ndim)=xc_tmp_recv(ix^d,1:ndim)*normconv(0)
2590 w_tec(1:nw+nwauxio)=wc_tmp_recv(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2591 write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2592 {end do\}
2593 end do
2594 case('tecplotCCmpi')
2595 ! in this option, we store the corner coordinates, and the cell center
2596 ! values of all variables. Due to this mix of corner/cell center, we must
2597 ! use BLOCK packaging, and thus we have enormous overhead by using
2598 ! calc_grid repeatedly to merely fill values of cell corner coordinates
2599 ! and cell center values per dimension, per variable
2600 if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
2601 if(nw+nwauxio==1)then
2602 ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
2603 ! and just set [ndim+1]
2604 if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2605 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
2606 'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2607 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2608 ndim+1,']=CELLCENTERED), ZONETYPE=', &
2609 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2610 else
2611 if(ndim+nw+nwauxio<10) then
2612 ! difference only in length of integer format specification for ndim+nw+nwauxio
2613 if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2614 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
2615 'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2616 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2617 ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2618 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2619 else
2620 if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2621 write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
2622 'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2623 ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2624 ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2625 {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2626 end if
2627 end if
2628
2629 do idim=1,ndim
2630 do morton_no=morton_start(ipe),morton_stop(ipe)
2631 itag=morton_no*idim
2632 call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2633 itag=igrid_recv*idim
2634 call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2635 if (level_recv/=level) cycle
2636 siz_ind=2*^nd
2637 itag=igrid_recv*idim
2638 call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2639 ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
2640 call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2641 ,icomm,intstatus(:,1),ierrmpi)
2642 call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2643 write(qunit,fmt="(100(e14.6))") xc_tmp_recv(ixrvc^s,idim)*normconv(0)
2644 end do
2645 end do
2646 do iw=1,nw+nwauxio
2647 do morton_no=morton_start(ipe),morton_stop(ipe)
2648 itag=morton_no*(ndim+iw)
2649 call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2650 itag=igrid_recv*(ndim+iw)
2651 call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2652 if (level_recv/=level) cycle
2653 siz_ind=2*^nd
2654 itag=igrid_recv*(ndim+iw)
2655 call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2656 ixrvccmin^d=ind_recv(^d);ixrvccmax^d=ind_recv(^nd+^d);
2657 call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2658 ,icomm,intstatus(:,1),ierrmpi)
2659 call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2660 write(qunit,fmt="(100(e14.6))") wcc_tmp_recv(ixrvcc^s,iw)*normconv(iw)
2661 end do
2662 end do
2663 case default
2664 call mpistop('no such tecplot type')
2665 end select
2666
2667 do morton_no=morton_start(ipe),morton_stop(ipe)
2668 itag=morton_no
2669 call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2670 itag=igrid_recv
2671 call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2672 if (level_recv/=level) cycle
2673 itag=igrid_recv
2674 call mpi_recv(igonlevel_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2675 call save_conntec(qunit,igrid_recv,igonlevel_recv)
2676 end do ! morton loop
2677 end do ! level loop
2678 end do ! ipe loop
2679 end if ! mype=0 if
2680 {^ifoned endif}
2681
2682 if (npe>1) then
2683 call mpi_barrier(icomm,ierrmpi)
2684 if(mype==0)deallocate(intstatus)
2685 end if
2686
2687 end subroutine tecplot_mpi
2688
2689 subroutine punstructuredvtkb_mpi(qunit)
2690 ! Write one pvtu and vtu files for each processor
2691 ! Otherwise like unstructuredvtk_mpi
2692 ! output for vtu format to paraview, binary version output
2693 ! uses calc_grid to compute nwauxio variables
2694 ! allows renormalizing using convert factors
2698
2699 integer, intent(in) :: qunit
2700
2701 double precision :: x_VTK(1:3)
2702 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
2703 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
2704 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
2705 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
2706 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
2707 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
2708 double precision :: normconv(0:nw+nwauxio)
2709 integer*8 :: offset
2710 integer :: igrid,iigrid,level,igonlevel,icel,ixC^L,ixCC^L,Morton_no
2711 integer :: NumGridsOnLevel(1:nlevelshi)
2712 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,VTK_type,ix^D
2713 integer:: recsep,k,iw,filenr
2714 integer:: length,lengthcc,offset_points,offset_cells, &
2715 length_coords,length_conn,length_offsets
2716 logical :: fileopen
2717 character:: buf
2718 character(len=6):: bufform
2719 character(len=80) :: pfilename
2720 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
2721 character(len=1024) :: outfilehead
2722
2723 ! Write pvtu-file:
2724 if (mype==0) then
2725 call write_pvtu(qunit)
2726 end if
2727 ! Now write the Source files:
2728 inquire(qunit,opened=fileopen)
2729 if(.not.fileopen)then
2730 ! generate filename
2731 filenr=snapshotnext-1
2732 if (autoconvert) filenr=snapshotnext
2733 ! Open the file for the header part
2734 write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename),filenr,"p",mype,".vtu"
2735 open(qunit,file=pfilename,status='unknown',form='formatted')
2736 end if
2737 ! generate xml header
2738 write(qunit,'(a)')'<?xml version="1.0"?>'
2739 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
2740 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
2741 write(qunit,'(a)')' <UnstructuredGrid>'
2742 write(qunit,'(a)')'<FieldData>'
2743 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
2744 'NumberOfTuples="1" format="ascii">'
2745 write(qunit,*) real(global_time*time_convert_factor)
2746 write(qunit,'(a)')'</DataArray>'
2747 write(qunit,'(a)')'</FieldData>'
2748 offset=0
2749 recsep=4
2750
2751 call getheadernames(wnamei,xandwnamei,outfilehead)
2752
2753 ! number of cells, number of corner points, per grid.
2754 nx^d=ixmhi^d-ixmlo^d+1;
2755 nxc^d=nx^d+1;
2756 nc={nx^d*}
2757 np={nxc^d*}
2758
2759 length=np*size_real
2760 lengthcc=nc*size_real
2761
2762 length_coords=3*length
2763 length_conn=2**^nd*size_int*nc
2764 length_offsets=nc*size_int
2765
2766 ! Note: using the w_write, writelevel, writespshift
2767 ! we can clip parts of the grid away, select variables, levels etc.
2768 do level=levmin,levmax
2769 if (writelevel(level)) then
2770 do morton_no=morton_start(mype),morton_stop(mype)
2771 igrid=sfc_to_igrid(morton_no)
2772 if (node(plevel_,igrid)/=level) cycle
2773 ! only output a grid when fully within clipped region selected
2774 ! by writespshift array
2775 if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
2776 *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
2777 <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
2778 select case(convert_type)
2779 case('pvtuBmpi')
2780 ! we write out every grid as one VTK PIECE
2781 write(qunit,'(a,i7,a,i7,a)') &
2782 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2783 write(qunit,'(a)')'<PointData>'
2784 do iw=1,nw
2785 if(.not.w_write(iw))cycle
2786
2787 write(qunit,'(a,a,a,i16,a)')&
2788 '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2789 '" format="appended" offset="',offset,'">'
2790 write(qunit,'(a)')'</DataArray>'
2791 offset=offset+length+size_int
2792 enddo
2793 do iw=nw+1,nw+nwauxio
2794
2795 write(qunit,'(a,a,a,i16,a)')&
2796 '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2797 '" format="appended" offset="',offset,'">'
2798 write(qunit,'(a)')'</DataArray>'
2799 offset=offset+length+size_int
2800 enddo
2801 write(qunit,'(a)')'</PointData>'
2802
2803 write(qunit,'(a)')'<Points>'
2804 write(qunit,'(a,i16,a)') &
2805 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
2806 ! write cell corner coordinates in a backward dimensional loop, always 3D output
2807 offset=offset+length_coords+size_int
2808 write(qunit,'(a)')'</Points>'
2809 case('pvtuBCCmpi')
2810 ! we write out every grid as one VTK PIECE
2811 write(qunit,'(a,i7,a,i7,a)') &
2812 '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2813 write(qunit,'(a)')'<CellData>'
2814 do iw=1,nw
2815 if(.not.w_write(iw))cycle
2816
2817 write(qunit,'(a,a,a,i16,a)')&
2818 '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2819 '" format="appended" offset="',offset,'">'
2820 write(qunit,'(a)')'</DataArray>'
2821 offset=offset+lengthcc+size_int
2822 enddo
2823 do iw=nw+1,nw+nwauxio
2824
2825 write(qunit,'(a,a,a,i16,a)')&
2826 '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2827 '" format="appended" offset="',offset,'">'
2828 write(qunit,'(a)')'</DataArray>'
2829 offset=offset+lengthcc+size_int
2830 enddo
2831 write(qunit,'(a)')'</CellData>'
2832 write(qunit,'(a)')'<Points>'
2833 write(qunit,'(a,i16,a)') &
2834 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
2835 ! write cell corner coordinates in a backward dimensional loop, always 3D output
2836 offset=offset+length_coords+size_int
2837 write(qunit,'(a)')'</Points>'
2838 end select
2839 write(qunit,'(a)')'<Cells>'
2840 ! connectivity part
2841 write(qunit,'(a,i16,a)')&
2842 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
2843 offset=offset+length_conn+size_int
2844 ! offsets data array
2845 write(qunit,'(a,i16,a)') &
2846 '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
2847 offset=offset+length_offsets+size_int
2848 ! VTK cell type data array
2849 write(qunit,'(a,i16,a)') &
2850 '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
2851 offset=offset+size_int+nc*size_int
2852 write(qunit,'(a)')'</Cells>'
2853 write(qunit,'(a)')'</Piece>'
2854 end if
2855 end do
2856 end if
2857 end do
2858
2859 write(qunit,'(a)')'</UnstructuredGrid>'
2860 write(qunit,'(a)')'<AppendedData encoding="raw">'
2861 close(qunit)
2862 ! next to make gfortran compiler happy, as it does not know
2863 ! form='binary' and produces error on compilation
2864 !bufform='binary'
2865 !open(qunit,file=pfilename,form=bufform,position='append')
2866 !This should in principle do also for gfortran (tested with gfortran 4.6.0 and Intel 11.1):
2867 open(qunit,file=pfilename,access='stream',form='unformatted',position='append')
2868 buf='_'
2869 write(qunit) trim(buf)
2870
2871 do level=levmin,levmax
2872 if (writelevel(level)) then
2873 do morton_no=morton_start(mype),morton_stop(mype)
2874 igrid=sfc_to_igrid(morton_no)
2875 if (node(plevel_,igrid)/=level) cycle
2876 ! only output a grid when fully within clipped region selected
2877 ! by writespshift array
2878 if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
2879 *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
2880 <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
2881 call calc_x(igrid,xc,xcc)
2882 call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2883 ixc^l,ixcc^l,.true.)
2884 do iw=1,nw
2885 if(.not.w_write(iw))cycle
2886 select case(convert_type)
2887 case('pvtuBmpi')
2888 write(qunit) length
2889 write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
2890 case('pvtuBCCmpi')
2891 write(qunit) lengthcc
2892 write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
2893 end select
2894 enddo
2895 do iw=nw+1,nw+nwauxio
2896 select case(convert_type)
2897 case('pvtuBmpi')
2898 write(qunit) length
2899 write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
2900 case('pvtuBCCmpi')
2901 write(qunit) lengthcc
2902 write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
2903 end select
2904 enddo
2905 write(qunit) length_coords
2906 {do ix^db=ixcmin^db,ixcmax^db \}
2907 x_vtk(1:3)=zero;
2908 x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
2909 do k=1,3
2910 write(qunit) real(x_vtk(k))
2911 end do
2912 {end do \}
2913 write(qunit) length_conn
2914 {do ix^db=1,nx^db\}
2915 {^ifoned write(qunit)ix1-1,ix1 \}
2916 {^iftwod
2917 write(qunit)(ix2-1)*nxc1+ix1-1, &
2918 (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
2919 \}
2920 {^ifthreed
2921 write(qunit)&
2922 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
2923 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
2924 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
2925 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
2926 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
2927 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
2928 ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
2929 ix3*nxc2*nxc1+ ix2*nxc1+ix1
2930 \}
2931 {end do\}
2932 write(qunit) length_offsets
2933 do icel=1,nc
2934 write(qunit) icel*(2**^nd)
2935 end do
2936 {^ifoned vtk_type=3 \}
2937 {^iftwod vtk_type=8 \}
2938 {^ifthreed vtk_type=11 \}
2939 write(qunit) size_int*nc
2940 do icel=1,nc
2941 write(qunit) vtk_type
2942 end do
2943 end if
2944 end do
2945 end if
2946 end do
2947
2948 close(qunit)
2949 open(qunit,file=pfilename,status='unknown',form='formatted',position='append')
2950 write(qunit,'(a)')'</AppendedData>'
2951 write(qunit,'(a)')'</VTKFile>'
2952 close(qunit)
2953
2954 end subroutine punstructuredvtkb_mpi
2955 {^iftwod
2956 ! subroutines to convert 2.5D data to 3D data
2957 subroutine unstructuredvtkb23(qunit)
2958 ! output for vtu format to paraview, binary version output
2959 ! not parallel, uses calc_grid to compute nwauxio variables
2961 use mod_physics
2963
2964 integer, intent(in) :: qunit
2965
2966 double precision :: x_VTK(1:3)
2967 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,3) :: xC_TMP
2968 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,& 3) :: xCC_TMP
2969 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,nw+nwauxio) :: wC_TMP
2970 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw& +nwauxio) :: wCC_TMP
2971 double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:nw& +nwauxio) :: w
2972 double precision :: normconv(0:nw+nwauxio)
2973 double precision :: zlength
2974 double precision ::d3grid,zlengsc,zgridsc
2975 integer*8 :: offset
2976 integer:: igrid,iigrid,level,igonlevel,icel,ixCmin1,ixCmin2,&
2977 ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,&
2978 ixCCmax2,ixCCmax3
2979 integer:: NumGridsOnLevel(1:nlevelshi)
2980 integer :: nx1,nx2,nx3,nxC1,nxC2,nxC3,nodesonlevel,elemsonlevel,nc,np,&
2981 VTK_type,ix1,ix2,ix3
2982 integer :: size_length,recsep,k,iw
2983 integer :: length,lengthcc,offset_points,offset_cells, length_coords,&
2984 length_conn,length_offsets
2985 integer :: i3grid,n3grid
2986 logical :: fileopen
2987 character:: buffer
2988 character(len=6):: bufform
2989 character(len=80):: filename
2990 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:3+nw+nwauxio)
2991 character(len=1024) :: outfilehead
2992
2993 if(npe>1)then
2994 if(mype==0) print *,'unstructuredvtkB23 not parallel, use vtumpi'
2995 call mpistop('npe>1, unstructuredvtkB23')
2996 end if
2997
2998 offset=0
2999 recsep=4
3000 size_length=4
3001 inquire(qunit,opened=fileopen)
3002 if(.not.fileopen)then
3003 ! generate filename
3004 write(filename,'(a,a,i4.4,a)') trim(base_filename),"3D",snapshotini,".vtu"
3005 ! Open the file for the header part
3006 open(qunit,file=filename,status='replace')
3007 endif
3008 call getheadernames(wnamei,xandwnamei,outfilehead)
3009 ! generate xml header
3010 write(qunit,'(a)')'<?xml version="1.0"?>'
3011 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
3012 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
3013 write(qunit,'(a)')'<UnstructuredGrid>'
3014 write(qunit,'(a)')'<FieldData>'
3015 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
3016 'NumberOfTuples="1" format="ascii">'
3017 write(qunit,'(f10.2)') real(global_time*time_convert_factor)
3018 write(qunit,'(a)')'</DataArray>'
3019 write(qunit,'(a)')'</FieldData>'
3020
3021 ! number of cells, number of corner points, per grid.
3022 nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
3023 nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
3024 nc=nx1*nx2*nx3
3025 np=nxc1*nxc2*nxc3
3026
3027 length=np*size_real
3028 lengthcc=nc*size_real
3029
3030 length_coords=3*length
3031 length_conn=2**3*size_int*nc
3032 length_offsets=nc*size_int
3033
3034 ! Note: using the w_write, writelevel, writespshift
3035 ! we can clip parts of the grid away, select variables, levels etc.
3036 zgridsc=2.d0
3037 zlengsc=2.d0*zgridsc
3038 zlength=zlengsc*(xprobmax1-xprobmin1)
3039 do level=levmin,levmax
3040 if (writelevel(level)) then
3041 do iigrid=1,igridstail; igrid=igrids(iigrid);
3042 if (node(plevel_,igrid)/=level) cycle
3043 block=>ps(igrid)
3044 ! only output a grid when fully within clipped region selected
3045 ! by writespshift array
3046 if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
3047 *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
3048 +(xprobmax2-xprobmin2)*writespshift(2,1))&
3049 .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
3050 *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
3051 -xprobmin2)*writespshift(2,2))) then
3052 d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
3053 n3grid=nint(zlength/d3grid)
3054 do i3grid=1,n3grid !subcycles
3055 select case(convert_type)
3056 case('vtuB23')
3057 ! we write out every grid as one VTK PIECE
3058 write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3059 '" NumberOfCells="',nc,'">'
3060 write(qunit,'(a)')'<PointData>'
3061 do iw=1,nw
3062 if(.not.w_write(iw))cycle
3063 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3064 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3065 write(qunit,'(a)')'</DataArray>'
3066 offset=offset+length+size_int
3067 enddo
3068 if(nwauxio>0)then
3069 do iw=nw+1,nw+nwauxio
3070 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3071 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3072 write(qunit,'(a)')'</DataArray>'
3073 offset=offset+length+size_int
3074 enddo
3075 endif
3076 write(qunit,'(a)')'</PointData>'
3077
3078 write(qunit,'(a)')'<Points>'
3079 write(qunit,'(a,i16,a)') &
3080 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3081 offset,'"/>'
3082 ! write cell corner coordinates in a backward dimensional loop, always 3D output
3083 offset=offset+length_coords+size_int
3084 write(qunit,'(a)')'</Points>'
3085 case('vtuBCC23')
3086 ! we write out every grid as one VTK PIECE
3087 write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3088 '" NumberOfCells="',nc,'">'
3089 write(qunit,'(a)')'<CellData>'
3090 do iw=1,nw
3091 if(.not.w_write(iw))cycle
3092 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3093 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3094 write(qunit,'(a)')'</DataArray>'
3095 offset=offset+lengthcc+size_int
3096 enddo
3097 if(nwauxio>0)then
3098 do iw=nw+1,nw+nwauxio
3099 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3100 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3101 write(qunit,'(a)')'</DataArray>'
3102 offset=offset+lengthcc+size_int
3103 enddo
3104 endif
3105 write(qunit,'(a)')'</CellData>'
3106 write(qunit,'(a)')'<Points>'
3107 write(qunit,'(a,i16,a)') &
3108 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3109 offset,'"/>'
3110 ! write cell corner coordinates in a backward dimensional loop, always 3D output
3111 offset=offset+length_coords+size_int
3112 write(qunit,'(a)')'</Points>'
3113 end select
3114 write(qunit,'(a)')'<Cells>'
3115 ! connectivity part
3116 write(qunit,'(a,i16,a)')&
3117 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
3118 offset,'"/>'
3119 offset=offset+length_conn+size_int
3120 ! offsets data array
3121 write(qunit,'(a,i16,a)') &
3122 '<DataArray type="Int32" Name="offsets" format="appended" offset="',&
3123 offset,'"/>'
3124 offset=offset+length_offsets+size_int
3125 ! VTK cell type data array
3126 write(qunit,'(a,i16,a)') &
3127 '<DataArray type="Int32" Name="types" format="appended" offset="',&
3128 offset,'"/>'
3129 offset=offset+size_length+nc*size_int
3130 write(qunit,'(a)')'</Cells>'
3131 write(qunit,'(a)')'</Piece>'
3132 end do !subcycles
3133 end if
3134 end do
3135 end if
3136 end do
3137
3138 write(qunit,'(a)')'</UnstructuredGrid>'
3139 write(qunit,'(a)')'<AppendedData encoding="raw">'
3140 close(qunit)
3141 open(qunit,file=filename,form='unformatted',access='stream',status='old',position='append')
3142 buffer='_'
3143 write(qunit) trim(buffer)
3144
3145 do level=levmin,levmax
3146 if (writelevel(level)) then
3147 do iigrid=1,igridstail; igrid=igrids(iigrid);
3148 if (node(plevel_,igrid)/=level) cycle
3149 block=>ps(igrid)
3150 ! only output a grid when fully within clipped region selected
3151 ! by writespshift array
3152 if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
3153 *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
3154 +(xprobmax2-xprobmin2)*writespshift(2,1))&
3155 .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
3156 *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
3157 -xprobmin2)*writespshift(2,2))) then
3158 d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
3159 n3grid=nint(zlength/d3grid)
3160 ! In case primitives to be saved: use primitive subroutine
3161 ! extra layer around mesh only needed when storing corner values and averaging
3162 if(saveprim) then
3163 call phys_to_primitive(ixglo1,ixglo2,ixghi1,ixghi2,&
3164 ixglo1,ixglo2,ixghi1,ixghi2,ps(igrid)%w,ps(igrid)%x)
3165 endif
3166 ! using array w so that new output auxiliaries can be calculated by the user
3167 ! extend 2D data to 3D insuring variables are independent on the third coordinate
3168 do ix3=ixglo1,ixghi1
3169 w(ixglo1:ixghi1,ixglo2:ixghi2,ix3,1:nw)=ps(igrid)%w(ixglo1:ixghi1,&
3170 ixglo2:ixghi2,1:nw)
3171 end do
3172 do i3grid=1,n3grid !subcycles
3173 call calc_grid23(qunit,igrid,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
3174 ixcmin1,ixcmin2,ixcmin3,ixcmax1,ixcmax2,ixcmax3,ixccmin1,ixccmin2,&
3175 ixccmin3,ixccmax1,ixccmax2,ixccmax3,.true.,i3grid,d3grid,w,zlength,zgridsc)
3176 do iw=1,nw
3177 if(.not.w_write(iw))cycle
3178 select case(convert_type)
3179 case('vtuB23')
3180 write(qunit) length
3181 write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3182 =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3183 case('vtuBCC23')
3184 write(qunit) lengthcc
3185 write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3186 =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3187 =ixccmin3,ixccmax3)
3188 end select
3189 enddo
3190 if(nwauxio>0)then
3191 do iw=nw+1,nw+nwauxio
3192 select case(convert_type)
3193 case('vtuB23')
3194 write(qunit) length
3195 write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3196 =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3197 case('vtuBCC23')
3198 write(qunit) lengthcc
3199 write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3200 =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3201 =ixccmin3,ixccmax3)
3202 end select
3203 end do
3204 end if
3205 write(qunit) length_coords
3206 do ix3=ixcmin3,ixcmax3
3207 do ix2=ixcmin2,ixcmax2
3208 do ix1=ixcmin1,ixcmax1
3209 x_vtk(1:3)=zero;
3210 x_vtk(1:3)=xc_tmp(ix1,ix2,ix3,1:3)*normconv(0);
3211 do k=1,3
3212 write(qunit) real(x_vtk(k))
3213 end do
3214 end do
3215 end do
3216 end do
3217 write(qunit) length_conn
3218 do ix3=1,nx3
3219 do ix2=1,nx2
3220 do ix1=1,nx1
3221 write(qunit)&
3222 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
3223 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3224 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3225 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
3226 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
3227 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3228 ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3229 ix3*nxc2*nxc1+ ix2*nxc1+ix1
3230 end do
3231 end do
3232 end do
3233 write(qunit) length_offsets
3234 do icel=1,nc
3235 write(qunit) icel*(2**3)
3236 end do
3237 vtk_type=11
3238 write(qunit) size_int*nc
3239 do icel=1,nc
3240 write(qunit) vtk_type
3241 end do
3242 end do !subcycles
3243 end if
3244 end do
3245 end if
3246 end do
3247
3248 close(qunit)
3249 open(qunit,file=filename,status='unknown',form='formatted',position='append')
3250
3251 write(qunit,'(a)')'</AppendedData>'
3252 write(qunit,'(a)')'</VTKFile>'
3253 close(qunit)
3254
3255 end subroutine unstructuredvtkb23
3256
3257 subroutine unstructuredvtkbsym23(qunit)
3258 ! output for vtu format to paraview, binary version output
3259 ! not parallel, uses calc_grid to compute nwauxio variables
3260 ! use this subroutine when the physical domain is symmetric/asymmetric about (0,y,z)
3261 ! plane, xprobmin1=0 and the computational domain is a half of the physical domain
3263 use mod_physics
3265
3266 integer, intent(in) :: qunit
3267
3268 double precision :: x_VTK(1:3)
3269 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,3) :: xC_TMP
3270 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,& 3) :: xCC_TMP
3271 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,nw+nwauxio) :: wC_TMP
3272 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw& +nwauxio) :: wCC_TMP
3273 double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:nw& +nwauxio) :: w
3274 double precision :: normconv(0:nw+nwauxio)
3275 double precision ::d3grid,zlengsc,zgridsc
3276 double precision :: zlength
3277 integer*8 :: offset
3278 integer:: igrid,iigrid,level,igonlevel,icel,ixCmin1,ixCmin2,&
3279 ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,&
3280 ixCCmax2,ixCCmax3
3281 integer:: NumGridsOnLevel(1:nlevelshi)
3282 integer :: nx1,nx2,nx3,nxC1,nxC2,nxC3,nodesonlevel,elemsonlevel,nc,np,&
3283 VTK_type,ix1,ix2,ix3
3284 integer :: size_length,recsep,k,iw
3285 integer :: length,lengthcc,offset_points,offset_cells, length_coords,&
3286 length_conn,length_offsets
3287 integer :: i3grid,n3grid
3288 logical :: fileopen
3289 character(len=80):: filename
3290 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:3+nw+nwauxio)
3291 character(len=1024) :: outfilehead
3292 character:: buffer
3293 character(len=6):: bufform
3294
3295 if(npe>1)then
3296 if(mype==0) print *,'unstructuredvtkBsym23 not parallel, use vtumpi'
3297 call mpistop('npe>1, unstructuredvtkBsym23')
3298 end if
3299
3300 offset=0
3301 recsep=4
3302 size_length=4
3303
3304 inquire(qunit,opened=fileopen)
3305 if(.not.fileopen)then
3306 ! generate filename
3307 write(filename,'(a,a,i4.4,a)') trim(base_filename),"3D",snapshotini,".vtu"
3308 ! Open the file for the header part
3309 open(qunit,file=filename,status='unknown')
3310 end if
3311
3312 call getheadernames(wnamei,xandwnamei,outfilehead)
3313 ! generate xml header
3314 write(qunit,'(a)')'<?xml version="1.0"?>'
3315 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
3316 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
3317 write(qunit,'(a)')'<UnstructuredGrid>'
3318 write(qunit,'(a)')'<FieldData>'
3319 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
3320 'NumberOfTuples="1" format="ascii">'
3321 write(qunit,'(f10.2)') real(global_time*time_convert_factor)
3322 write(qunit,'(a)')'</DataArray>'
3323 write(qunit,'(a)')'</FieldData>'
3324
3325 ! number of cells, number of corner points, per grid.
3326 nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
3327 nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
3328 nc=nx1*nx2*nx3
3329 np=nxc1*nxc2*nxc3
3330
3331 length=np*size_real
3332 lengthcc=nc*size_real
3333
3334 length_coords=3*length
3335 length_conn=2**3*size_int*nc
3336 length_offsets=nc*size_int
3337
3338 ! Note: using the w_write, writelevel, writespshift
3339 ! we can clip parts of the grid away, select variables, levels etc.
3340 zlengsc=4.d0
3341 zgridsc=2.d0
3342 zlength=zlengsc*(xprobmax1-xprobmin1)
3343 do level=levmin,levmax
3344 if (writelevel(level)) then
3345 do iigrid=1,igridstail; igrid=igrids(iigrid);
3346 if (node(plevel_,igrid)/=level) cycle
3347 block=>ps(igrid)
3348 ! only output a grid when fully within clipped region selected
3349 ! by writespshift array
3350 if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
3351 *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
3352 +(xprobmax2-xprobmin2)*writespshift(2,1))&
3353 .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
3354 *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
3355 -xprobmin2)*writespshift(2,2))) then
3356 d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
3357 n3grid=nint(zlength/d3grid)
3358 do i3grid=1,n3grid !subcycles
3359 !! original domain ----------------------------------start
3360 select case(convert_type)
3361 case('vtuBsym23')
3362 ! we write out every grid as one VTK PIECE
3363 write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3364 '" NumberOfCells="',nc,'">'
3365 write(qunit,'(a)')'<PointData>'
3366 do iw=1,nw
3367 if(.not.w_write(iw))cycle
3368 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3369 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3370 write(qunit,'(a)')'</DataArray>'
3371 offset=offset+length+size_length
3372 enddo
3373 if(nwauxio>0)then
3374 do iw=nw+1,nw+nwauxio
3375 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3376 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3377 write(qunit,'(a)')'</DataArray>'
3378 offset=offset+length+size_length
3379 enddo
3380 endif
3381 write(qunit,'(a)')'</PointData>'
3382 write(qunit,'(a)')'<Points>'
3383 write(qunit,'(a,i16,a)') &
3384 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3385 offset,'"/>'
3386 ! write cell corner coordinates in a backward dimensional loop, always 3D output
3387 offset=offset+length_coords+size_length
3388 write(qunit,'(a)')'</Points>'
3389 case('vtuBCCsym23')
3390 ! we write out every grid as one VTK PIECE
3391 write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3392 '" NumberOfCells="',nc,'">'
3393 write(qunit,'(a)')'<CellData>'
3394 do iw=1,nw
3395 if(.not.w_write(iw))cycle
3396 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3397 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3398 write(qunit,'(a)')'</DataArray>'
3399 offset=offset+lengthcc+size_length
3400 enddo
3401 if(nwauxio>0)then
3402 do iw=nw+1,nw+nwauxio
3403 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3404 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3405 write(qunit,'(a)')'</DataArray>'
3406 offset=offset+lengthcc+size_length
3407 enddo
3408 endif
3409 write(qunit,'(a)')'</CellData>'
3410
3411 write(qunit,'(a)')'<Points>'
3412 write(qunit,'(a,i16,a)') &
3413 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3414 offset,'"/>'
3415 ! write cell corner coordinates in a backward dimensional loop, always 3D output
3416 offset=offset+length_coords+size_length
3417 write(qunit,'(a)')'</Points>'
3418 end select
3419 write(qunit,'(a)')'<Cells>'
3420 ! connectivity part
3421 write(qunit,'(a,i16,a)')&
3422 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
3423 offset,'"/>'
3424 offset=offset+length_conn+size_length
3425 ! offsets data array
3426 write(qunit,'(a,i16,a)') &
3427 '<DataArray type="Int32" Name="offsets" format="appended" offset="',&
3428 offset,'"/>'
3429 offset=offset+length_offsets+size_length
3430 ! VTK cell type data array
3431 write(qunit,'(a,i16,a)') &
3432 '<DataArray type="Int32" Name="types" format="appended" offset="',&
3433 offset,'"/>'
3434 offset=offset+size_length+nc*size_int
3435 write(qunit,'(a)')'</Cells>'
3436 write(qunit,'(a)')'</Piece>'
3437 !! original domain ----------------------------------end
3438 !! symetric/asymetric mirror domain -----------------start
3439 select case(convert_type)
3440 case('vtuBsym23')
3441 ! we write out every grid as one VTK PIECE
3442 write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3443 '" NumberOfCells="',nc,'">'
3444 write(qunit,'(a)')'<PointData>'
3445 do iw=1,nw
3446 if(.not.w_write(iw))cycle
3447 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3448 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3449 write(qunit,'(a)')'</DataArray>'
3450 offset=offset+length+size_length
3451 enddo
3452 if(nwauxio>0)then
3453 do iw=nw+1,nw+nwauxio
3454 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3455 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3456 write(qunit,'(a)')'</DataArray>'
3457 offset=offset+length+size_length
3458 enddo
3459 endif
3460 write(qunit,'(a)')'</PointData>'
3461 write(qunit,'(a)')'<Points>'
3462 write(qunit,'(a,i16,a)') &
3463 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3464 offset,'"/>'
3465 ! write cell corner coordinates in a backward dimensional loop, always 3D output
3466 offset=offset+length_coords+size_length
3467 write(qunit,'(a)')'</Points>'
3468 case('vtuBCCsym23')
3469 ! we write out every grid as one VTK PIECE
3470 write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3471 '" NumberOfCells="',nc,'">'
3472 write(qunit,'(a)')'<CellData>'
3473 do iw=1,nw
3474 if(.not.w_write(iw))cycle
3475 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3476 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3477 write(qunit,'(a)')'</DataArray>'
3478 offset=offset+lengthcc+size_length
3479 enddo
3480 if(nwauxio>0)then
3481 do iw=nw+1,nw+nwauxio
3482 write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3483 trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3484 write(qunit,'(a)')'</DataArray>'
3485 offset=offset+lengthcc+size_length
3486 enddo
3487 endif
3488 write(qunit,'(a)')'</CellData>'
3489 write(qunit,'(a)')'<Points>'
3490 write(qunit,'(a,i16,a)') &
3491 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3492 offset,'"/>'
3493 ! write cell corner coordinates in a backward dimensional loop, always 3D output
3494 offset=offset+length_coords+size_length
3495 write(qunit,'(a)')'</Points>'
3496 end select
3497 write(qunit,'(a)')'<Cells>'
3498 ! connectivity part
3499 write(qunit,'(a,i16,a)')&
3500 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
3501 offset,'"/>'
3502 offset=offset+length_conn+size_length
3503 ! offsets data array
3504 write(qunit,'(a,i16,a)') &
3505 '<DataArray type="Int32" Name="offsets" format="appended" offset="',&
3506 offset,'"/>'
3507 offset=offset+length_offsets+size_length
3508 ! VTK cell type data array
3509 write(qunit,'(a,i16,a)') &
3510 '<DataArray type="Int32" Name="types" format="appended" offset="',&
3511 offset,'"/>'
3512 offset=offset+size_length+nc*size_int
3513 write(qunit,'(a)')'</Cells>'
3514 write(qunit,'(a)')'</Piece>'
3515 !! symetric/asymetric mirror domain -----------------end
3516 end do !subcycles
3517 end if
3518 end do
3519 end if
3520 end do
3521
3522 write(qunit,'(a)')'</UnstructuredGrid>'
3523 write(qunit,'(a)')'<AppendedData encoding="raw">'
3524 close(qunit)
3525 open(qunit,file=filename,form='unformatted',access='stream',status='old',position='append')
3526 buffer='_'
3527 write(qunit) trim(buffer)
3528 do level=levmin,levmax
3529 if (writelevel(level)) then
3530 do iigrid=1,igridstail; igrid=igrids(iigrid);
3531 if (node(plevel_,igrid)/=level) cycle
3532 block=>ps(igrid)
3533 ! only output a grid when fully within clipped region selected
3534 ! by writespshift array
3535 if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
3536 *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
3537 +(xprobmax2-xprobmin2)*writespshift(2,1))&
3538 .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
3539 *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
3540 -xprobmin2)*writespshift(2,2))) then
3541 d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
3542 n3grid=nint(zlength/d3grid)
3543 ! In case primitives to be saved: use primitive subroutine
3544 ! extra layer around mesh only needed when storing corner values and averaging
3545 if(saveprim) then
3546 call phys_to_primitive(ixglo1,ixglo2,ixghi1,ixghi2,&
3547 ixglo1,ixglo2,ixghi1,ixghi2,ps(igrid)%w,ps(igrid)%x)
3548 endif
3549 ! using array w so that new output auxiliaries can be calculated by the user
3550 ! extend 2D data to 3D insuring variables are independent on the third coordinate
3551 do ix3=ixglo1,ixghi1
3552 w(ixglo1:ixghi1,ixglo2:ixghi2,ix3,1:nw)=ps(igrid)%w(ixglo1:ixghi1,&
3553 ixglo2:ixghi2,1:nw)
3554 end do
3555 do i3grid=1,n3grid !subcycles
3556 call calc_grid23(qunit,igrid,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
3557 ixcmin1,ixcmin2,ixcmin3,ixcmax1,ixcmax2,ixcmax3,ixccmin1,ixccmin2,&
3558 ixccmin3,ixccmax1,ixccmax2,ixccmax3,.true.,i3grid,d3grid,w,zlength,zgridsc)
3559 !! original domain ----------------------------------start
3560 do iw=1,nw
3561 if(.not.w_write(iw))cycle
3562 select case(convert_type)
3563 case('vtuBsym23')
3564 write(qunit) length
3565 write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3566 =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3567 case('vtuBCCsym23')
3568 write(qunit) lengthcc
3569 write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3570 =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3571 =ixccmin3,ixccmax3)
3572 end select
3573 enddo
3574 if(nwauxio>0)then
3575 do iw=nw+1,nw+nwauxio
3576 select case(convert_type)
3577 case('vtuBsym23')
3578 write(qunit) length
3579 write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3580 =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3581 case('vtuBCCsym23')
3582 write(qunit) lengthcc
3583 write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3584 =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3585 =ixccmin3,ixccmax3)
3586 end select
3587 enddo
3588 endif
3589 write(qunit) length_coords
3590 do ix3=ixcmin3,ixcmax3
3591 do ix2=ixcmin2,ixcmax2
3592 do ix1=ixcmin1,ixcmax1
3593 x_vtk(1:3)=zero;
3594 x_vtk(1:3)=xc_tmp(ix1,ix2,ix3,1:3)*normconv(0);
3595 do k=1,3
3596 write(qunit) real(x_vtk(k))
3597 end do
3598 end do
3599 end do
3600 end do
3601 write(qunit) length_conn
3602 do ix3=1,nx3
3603 do ix2=1,nx2
3604 do ix1=1,nx1
3605 write(qunit)&
3606 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
3607 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3608 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3609 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
3610 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
3611 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3612 ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3613 ix3*nxc2*nxc1+ ix2*nxc1+ix1
3614 end do
3615 end do
3616 end do
3617 write(qunit) length_offsets
3618 do icel=1,nc
3619 write(qunit) icel*(2**3)
3620 end do
3621 vtk_type=11
3622 write(qunit) size_int*nc
3623 do icel=1,nc
3624 write(qunit) vtk_type
3625 end do
3626 !! original domain ----------------------------------end
3627 !! symetric/asymetric mirror domain -----------------start
3628 do iw=1,nw
3629 if(.not.w_write(iw))cycle
3630 if(iw==2 .or. iw==4 .or. iw==7) then
3631 wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,iw)=&
3632 -wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,iw)
3633 wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,iw)=&
3634 -wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,iw)
3635 end if
3636 select case(convert_type)
3637 case('vtuBsym23')
3638 write(qunit) length
3639 write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3640 =ixcmax1,ixcmin1,-1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3641 case('vtuBCCsym23')
3642 write(qunit) lengthcc
3643 write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3644 =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3645 =ixccmin3,ixccmax3)
3646 end select
3647 enddo
3648 if(nwauxio>0)then
3649 do iw=nw+1,nw+nwauxio
3650 select case(convert_type)
3651 case('vtuBsym23')
3652 write(qunit) length
3653 write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3654 =ixcmax1,ixcmin1,-1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3655 case('vtuBCCsym23')
3656 write(qunit) lengthcc
3657 write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3658 =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3659 =ixccmin3,ixccmax3)
3660 end select
3661 end do
3662 end if
3663 write(qunit) length_coords
3664 do ix3=ixcmin3,ixcmax3
3665 do ix2=ixcmin2,ixcmax2
3666 do ix1=ixcmax1,ixcmin1,-1
3667 x_vtk(1:3)=zero;
3668 x_vtk(1:3)=xc_tmp(ix1,ix2,ix3,1:3)*normconv(0);
3669 x_vtk(1)=-x_vtk(1)
3670 do k=1,3
3671 write(qunit) real(x_vtk(k))
3672 end do
3673 end do
3674 end do
3675 end do
3676 write(qunit) length_conn
3677 do ix3=1,nx3
3678 do ix2=1,nx2
3679 do ix1=nx1,1,-1
3680 write(qunit)&
3681 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
3682 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3683 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3684 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
3685 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
3686 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3687 ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3688 ix3*nxc2*nxc1+ ix2*nxc1+ix1
3689 end do
3690 end do
3691 end do
3692 write(qunit) length_offsets
3693 do icel=1,nc
3694 write(qunit) icel*(2**3)
3695 end do
3696 vtk_type=11
3697 write(qunit) size_int*nc
3698 do icel=1,nc
3699 write(qunit) vtk_type
3700 end do
3701 !! symetric/asymetric mirror domain -----------------end
3702 end do !subcycles
3703 end if
3704 end do
3705 end if
3706 end do
3707 close(qunit)
3708 open(qunit,file=filename,status='unknown',form='formatted',position='append')
3709 write(qunit,'(a)')'</AppendedData>'
3710 write(qunit,'(a)')'</VTKFile>'
3711 close(qunit)
3712
3713 end subroutine unstructuredvtkbsym23
3714
3715 subroutine calc_grid23(qunit,igrid,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP,normconv,&
3716 ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,&
3717 ixCCmax1,ixCCmax2,ixCCmax3,first,i3grid,d3grid,w,zlength,zgridsc)
3718 ! this subroutine computes both corner as well as cell-centered values
3719 ! it handles how we do the center to corner averaging, as well as
3720 ! whether we switch to cartesian or want primitive or conservative output,
3721 ! handling the addition of B0 in B0+B1 cases, ...
3722 ! the normconv is passed on to specialvar_output for extending with
3723 ! possible normalization values for the nw+1:nw+nwauxio entries
3725 integer, intent(in) :: qunit, igrid,i3grid
3726 logical, intent(in) :: first
3727
3728 double precision :: dx1,dx2,dx3,d3grid,zlength,zgridsc
3729 double precision :: ldw(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1),&
3730 dwC(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1)
3731 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,3) :: xC
3732 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,& 3) :: xCC
3733 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,nw+nwauxio) :: wC
3734 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw& +nwauxio) :: wCC
3735 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,3) :: xC_TMP
3736 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,& 3) :: xCC_TMP
3737 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,nw+nwauxio) :: wC_TMP
3738 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw& +nwauxio) :: wCC_TMP
3739 double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:nw& +nwauxio) :: w
3740 double precision,dimension(0:nw+nwauxio) :: normconv
3741 integer :: nx1,nx2,nx3, nxC1,nxC2,nxC3, ix1,ix2,ix3, ix, iw, level, idir
3742 integer :: ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,&
3743 ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3,nxCC1,nxCC2,nxCC3
3744 integer :: idims,jxCmin1,jxCmin2,jxCmin3,jxCmax1,jxCmax2,jxCmax3
3745 logical, save :: subfirst=.true.
3746
3747 ! following only for allowing compiler to go through with debug on
3748 nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
3749 level=node(plevel_,igrid)
3750 dx1=dx(1,level);dx2=dx(2,level);dx3=zgridsc*dx(1,level);
3751 ! for normalization within the code
3752 if(saveprim) then
3753 normconv(0) = length_convert_factor
3754 normconv(1:nw) = w_convert_factor
3755 else
3756 normconv(0)=length_convert_factor
3757 ! assuming density
3758 normconv(1)=w_convert_factor(1)
3759 ! assuming momentum=density*velocity
3760 if (nw>=2) normconv(2:2+3)=w_convert_factor(1)*w_convert_factor(2:2+3)
3761 ! assuming energy/pressure and magnetic field
3762 if (nw>=2+3) normconv(2+3:nw)=w_convert_factor(2+3:nw)
3763 end if
3764 ! coordinates of cell centers
3765 nxcc1=nx1;nxcc2=nx2;nxcc3=nx3;
3766 ixccmin1=ixmlo1;ixccmin2=ixmlo2;ixccmin3=ixmlo1; ixccmax1=ixmhi1
3767 ixccmax2=ixmhi2;ixccmax3=ixmhi1;
3768 do ix=ixccmin1,ixccmax1
3769 xcc(ix,ixccmin2:ixccmax2,ixccmin3:ixccmax3,1)=rnode(rpxmin1_,igrid)&
3770 +(dble(ix-ixccmin1)+half)*dx1
3771 end do
3772 do ix=ixccmin2,ixccmax2
3773 xcc(ixccmin1:ixccmax1,ix,ixccmin3:ixccmax3,2)=rnode(rpxmin2_,igrid)&
3774 +(dble(ix-ixccmin2)+half)*dx2
3775 end do
3776 do ix=ixccmin3,ixccmax3
3777 xcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ix,3)=-zlength/two+&
3778 dble(i3grid-1)*d3grid+(dble(ix-ixccmin3)+half)*dx3
3779 end do
3780
3781 ! coordinates of cell corners
3782 nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
3783 ixcmin1=ixmlo1-1;ixcmin2=ixmlo2-1;ixcmin3=ixmlo1-1; ixcmax1=ixmhi1
3784 ixcmax2=ixmhi2;ixcmax3=ixmhi1;
3785 do ix=ixcmin1,ixcmax1
3786 xc(ix,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1)=rnode(rpxmin1_,igrid)&
3787 +dble(ix-ixcmin1)*dx1
3788 end do
3789 do ix=ixcmin2,ixcmax2
3790 xc(ixcmin1:ixcmax1,ix,ixcmin3:ixcmax3,2)=rnode(rpxmin2_,igrid)&
3791 +dble(ix-ixcmin2)*dx2
3792 end do
3793 do ix=ixcmin3,ixcmax3
3794 xc(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ix,3)=-zlength/two+&
3795 dble(i3grid-1)*d3grid+dble(ix-ixcmin3)*dx3
3796 end do
3797
3798 if (nwextra>0) then
3799 ! here we actually fill the ghost layers for the nwextra variables using
3800 ! continuous extrapolation (as these values do not exist normally in ghost
3801 ! cells)
3802 do idims=1,3
3803 select case(idims)
3804 case(1)
3805 jxcmin1=ixghi1+1-nghostcells;jxcmin2=ixglo2;jxcmin3=ixglo1;
3806 jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixghi1;
3807 do ix1=jxcmin1,jxcmax1
3808 w(ix1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw) = w(jxcmin1&
3809 -1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw)
3810 end do
3811 jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixglo1;
3812 jxcmax1=ixglo1-1+nghostcells;jxcmax2=ixghi2;jxcmax3=ixghi1;
3813 do ix1=jxcmin1,jxcmax1
3814 w(ix1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw) = w(jxcmax1&
3815 +1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw)
3816 end do
3817 case(2)
3818 jxcmin1=ixglo1;jxcmin2=ixghi2+1-nghostcells;jxcmin3=ixglo1;
3819 jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixghi1;
3820 do ix2=jxcmin2,jxcmax2
3821 w(jxcmin1:jxcmax1,ix2,jxcmin3:jxcmax3,nw-nwextra+1:nw) &
3822 = w(jxcmin1:jxcmax1,jxcmin2-1,jxcmin3:jxcmax3,nw-nwextra+1:nw)
3823 end do
3824 jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixglo1;
3825 jxcmax1=ixghi1;jxcmax2=ixglo2-1+nghostcells;jxcmax3=ixghi1;
3826 do ix2=jxcmin2,jxcmax2
3827 w(jxcmin1:jxcmax1,ix2,jxcmin3:jxcmax3,nw-nwextra+1:nw) &
3828 = w(jxcmin1:jxcmax1,jxcmax2+1,jxcmin3:jxcmax3,nw-nwextra+1:nw)
3829 end do
3830 case(3)
3831 jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixghi1+1-nghostcells;
3832 jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixghi1;
3833 do ix3=jxcmin3,jxcmax3
3834 w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,ix3,nw-nwextra+1:nw) &
3835 = w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,jxcmin3-1,nw-nwextra+1:nw)
3836 end do
3837 jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixglo1;
3838 jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixglo1-1+nghostcells;
3839 do ix3=jxcmin3,jxcmax3
3840 w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,ix3,nw-nwextra+1:nw) &
3841 = w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,jxcmax3+1,nw-nwextra+1:nw)
3842 end do
3843 end select
3844 end do
3845 end if
3846 ! next lines needed when specialvar_output uses gradients
3847 ! and later on when dwlimiter2 is used
3848 if(nwauxio>0)then
3849 ! auxiliary io variables can be computed and added by user
3850 ! next few lines ensure correct usage of routines like divvector etc
3851 dxlevel(1)=rnode(rpdx1_,igrid);dxlevel(2)=rnode(rpdx2_,igrid)
3852 ! default (no) normalization for auxiliary variables
3853 normconv(nw+1:nw+nwauxio)=one
3854 ! maybe need for restriction to ixG^LL^LSUB1
3855 call specialvar_output23(ixglo1,ixglo2,ixglo1,ixghi1,ixghi2,ixghi1,ixglo1&
3856 +1,ixglo2+1,ixglo1+1,ixghi1-1,ixghi2-1,ixghi1-1,w,xcc,normconv)
3857 endif
3858 ! compute the cell-center values for w first
3859 !===========================================
3860 ! cell center values obtained from mere copy, while B0+B1 split handled here
3861 wcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,:)=w(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,:)
3862 if(b0field) then
3863 do ix3=ixccmin3,ixccmax3
3864 do ix2=ixccmin2,ixccmax2
3865 do ix1=ixccmin1,ixccmax1
3866 wcc(ix1,ix2,ix3,iw_mag(:))=wcc(ix1,ix2,ix3,iw_mag(:))+ps(igrid)%B0(ix1,ix2,&
3867 :,0)
3868 end do
3869 end do
3870 end do
3871 end if
3872 if(.not.saveprim .and. b0field .and. iw_e>0) then
3873 do ix3=ixccmin3,ixccmax3
3874 do ix2=ixccmin2,ixccmax2
3875 do ix1=ixccmin1,ixccmax1
3876 wcc(ix1,ix2,ix3,iw_e)=w(ix1,ix2,ix3,iw_e) +half*sum(ps(igrid)%B0(ix1,&
3877 ix2,:,0)**2 ) + sum(w(ix1,ix2,ix3,&
3878 iw_mag(:))*ps(igrid)%B0(ix1,ix2,:,0))
3879 end do
3880 end do
3881 end do
3882 end if
3883 ! compute the corner values for w now by averaging
3884 !=================================================
3885 if(slab_uniform)then
3886 ! for slab symmetry: no geometrical info required
3887 do iw=1,nw+nwauxio
3888 if (b0field.and.iw>iw_mag(1)-1.and.iw<=iw_mag(ndir)) then
3889 idir=iw-iw_mag(1)+1
3890 do ix3=ixcmin3,ixcmax3
3891 do ix2=ixcmin2,ixcmax2
3892 do ix1=ixcmin1,ixcmax1
3893 wc(ix1,ix2,ix3,iw)=sum(w(ix1:ix1+1,ix2:ix2+1,ix3,iw) &
3894 +ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
3895 ,idir,0))/dble(2**3)+&
3896 sum(w(ix1:ix1+1,ix2:ix2+1,ix3+1,iw) &
3897 +ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
3898 ,idir,0))/dble(2**3)
3899 end do
3900 end do
3901 end do
3902 else
3903 do ix3=ixcmin3,ixcmax3
3904 do ix2=ixcmin2,ixcmax2
3905 do ix1=ixcmin1,ixcmax1
3906 wc(ix1,ix2,ix3,iw)=sum(w(ix1:ix1+1,ix2:ix2+1,ix3:ix3&
3907 +1,iw))/dble(2**3)
3908 end do
3909 end do
3910 end do
3911 end if
3912 end do
3913 if(.not.saveprim .and. b0field .and. iw_e>0) then
3914 do ix3=ixcmin3,ixcmax3
3915 do ix2=ixcmin2,ixcmax2
3916 do ix1=ixcmin1,ixcmax1
3917 wc(ix1,ix2,ix3,iw_e)=sum( w(ix1:ix1+1,ix2:ix2+1,ix3,iw_e) &
3918 +half*sum(ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
3919 ,:,0)**2,dim=ndim+1) + sum( w(ix1:ix1+1,ix2:ix2+1,ix3&
3920 ,iw_mag(:))*ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
3921 ,:,0),dim=ndim+1) ) /dble(2**3)+&
3922 sum( w(ix1:ix1+1,ix2:ix2+1,ix3+1,iw_e) &
3923 +half*sum(ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
3924 ,:,0)**2,dim=ndim+1) + sum( w(ix1:ix1+1,ix2:ix2+1,ix3&
3925 +1,iw_mag(:))*ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
3926 ,:,0),dim=ndim+1) ) /dble(2**3)
3927 end do
3928 end do
3929 end do
3930 end if
3931 end if
3932 ! keep the coordinate and vector components
3933 xc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:3) &
3934 = xc(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:3)
3935 wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:nw&
3936 +nwauxio) = wc(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:nw&
3937 +nwauxio)
3938 xcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,&
3939 1:3) = xcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,&
3940 ixccmin3:ixccmax3,1:3)
3941 wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:nw&
3942 +nwauxio) = wcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,&
3943 1:nw+nwauxio)
3944 end subroutine calc_grid23
3945
3946 subroutine save_connvtk23(qunit,igrid)
3947 ! this saves the basic line, pixel and voxel connectivity,
3948 ! as used by VTK file outputs for unstructured grid
3950
3951 integer, intent(in) :: qunit, igrid
3952
3953 integer :: nx1,nx2,nx3, nxC1,nxC2,nxC3, ix1,ix2,ix3
3954
3955 nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
3956 nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
3957 do ix3=1,nx3
3958 do ix2=1,nx2
3959 do ix1=1,nx1
3960 write(qunit,'(8(i7,1x))')&
3961 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
3962 (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3963 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3964 (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
3965 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
3966 ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3967 ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3968 ix3*nxc2*nxc1+ ix2*nxc1+ix1
3969
3970 end do
3971 end do
3972 end do
3973 end subroutine save_connvtk23
3974
3975 subroutine specialvar_output23(ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
3976 ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,w,x,normconv)
3977 ! this subroutine can be used in convert, to add auxiliary variables to the
3978 ! converted output file, for further analysis using tecplot, paraview, ....
3979 ! these auxiliary values need to be stored in the nw+1:nw+nwauxio slots
3980 ! the array normconv can be filled in the (nw+1:nw+nwauxio) range with
3981 ! corresponding normalization values (default value 1)
3983
3984 integer, intent(in) :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
3985 ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
3986 double precision, intent(in) :: x(ixImin1:ixImax1,ixImin2:ixImax2,&
3987 ixImin3:ixImax3,1:3)
3988 double precision :: w(ixImin1:ixImax1,ixImin2:ixImax2,&
3989 ixImin3:ixImax3,nw+nwauxio)
3990 double precision :: normconv(0:nw+nwauxio)
3991
3992 double precision :: qvec(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:ndir),&
3993 curlvec(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:ndir)
3994 integer :: idirmin
3995
3996 ! output Te
3997 !if(saveprim)then
3998 ! w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+1)=w(ixOmin1:ixOmax1,&
3999 ! ixOmin2:ixOmax2,ixOmin3:ixOmax3,iw_e)/w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
4000 ! ixOmin3:ixOmax3,iw_rho)
4001 !endif
4002 !!! store current
4003 ! qvec(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3,1)=w(ixImin1:ixImax1,&
4004 ! ixImin2:ixImax2,ixImin3:ixImax3,mag(1))
4005 ! qvec(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3,2)=w(ixImin1:ixImax1,&
4006 ! ixImin2:ixImax2,ixImin3:ixImax3,mag(2))
4007 ! qvec(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3,3)=w(ixImin1:ixImax1,&
4008 ! ixImin2:ixImax2,ixImin3:ixImax3,mag(3));
4009 !call curlvector3D(qvec,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,ixOmin1,&
4010 ! ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,curlvec,idirmin,1,ndir)
4011 !w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+2)=curlvec&
4012 ! (ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,1)
4013 !w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+3)=curlvec&
4014 ! (ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,2)
4015 !w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+4)=curlvec&
4016 ! (ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,3);
4017 end subroutine specialvar_output23
4018 \}
4019
4020end module mod_convert_files
4021
subroutine, public alloc_state_output(igrid, s, ixgl)
allocate memory to physical state of igrid node for output with auxio
Handles computations for coordinates and variables in output.
subroutine getheadernames(wnamei, xandwnamei, outfilehead)
get all variables names
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.
subroutine calc_x(igrid, xc, xcc)
computes cell corner (xC) and cell center (xCC) coordinates
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
subroutine unstructuredvtkb64(qunit)
subroutine punstructuredvtk_mpi(qunit)
subroutine write_vtk(qunit, ixil, ixcl, ixccl, igrid, nc, np, nxd, nxcd, normconv, wnamei, xc, xcc, wc, wcc)
subroutine tecplot_mpi(qunit)
subroutine oneblock(qunit)
subroutine write_pvtu(qunit)
subroutine imagedatavtk_mpi(qunit)
subroutine save_connvtk(qunit, igrid)
integer function nodenumbertec2d(i1, i2, nx1, nx2, ig, igrid)
subroutine calc_grid23(qunit, igrid, xc_tmp, xcc_tmp, wc_tmp, wcc_tmp, normconv, ixcmin1, ixcmin2, ixcmin3, ixcmax1, ixcmax2, ixcmax3, ixccmin1, ixccmin2, ixccmin3, ixccmax1, ixccmax2, ixccmax3, first, i3grid, d3grid, w, zlength, zgridsc)
subroutine tecplot(qunit)
subroutine generate_plotfile
integer function nodenumbertec1d(i1, nx1, ig, igrid)
subroutine unstructuredvtk(qunit)
subroutine unstructuredvtkb23(qunit)
subroutine unstructuredvtkb(qunit)
subroutine write_vti(qunit, ixil, ixcl, ixccl, igd, nxd, normconv, wnamei, wc, wcc)
integer function nodenumbertec3d(i1, i2, i3, nx1, nx2, nx3, ig, igrid)
subroutine save_conntec(qunit, igrid, igonlevel)
subroutine save_connvtk23(qunit, igrid)
subroutine onegrid(qunit)
subroutine unstructuredvtk_mpi(qunit)
subroutine punstructuredvtkb_mpi(qunit)
subroutine unstructuredvtkbsym23(qunit)
subroutine specialvar_output23(iximin1, iximin2, iximin3, iximax1, iximax2, iximax3, ixomin1, ixomin2, ixomin3, ixomax1, ixomax2, ixomax3, w, x, normconv)
subroutine convert_all()
Definition mod_convert.t:91
Module with basic grid data structures.
Definition mod_forest.t:2
integer, dimension(:), allocatable, save sfc_to_igrid
Go from a Morton number to an igrid index (for a single processor)
Definition mod_forest.t:53
integer, dimension(:), allocatable, save morton_start
First Morton number per processor.
Definition mod_forest.t:62
integer, dimension(:), allocatable, save morton_stop
Last Morton number per processor.
Definition mod_forest.t:65
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
update ghost cells of all blocks including physical boundaries
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision, dimension(:), allocatable w_convert_factor
Conversion factors the primitive variables.
type(state), pointer block
Block pointer for using one block and its previous state.
logical nocartesian
IO switches for conversion.
double precision global_time
The global simulation time.
integer type_block_xc_io
MPI type for IO: cell corner (xc) or cell center (xcc) coordinates.
integer snapshotini
Resume from the snapshot with this index.
logical saveprim
If true, convert from conservative to primitive variables in output.
character(len=std_len) convert_type
Which format to use when converting.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision time_convert_factor
Conversion factor for time unit.
integer icomm
The MPI communicator.
integer, dimension(:), allocatable ng
number of grid blocks in domain per dimension, in array over levels
integer mype
The rank of the current MPI task.
integer type_block_io
MPI type for IO: block excluding ghost cells.
double precision length_convert_factor
Conversion factor for length unit.
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
integer type_block_wc_io
MPI type for IO: cell corner (wc) or cell center (wcc) variables.
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
integer npe
The number of MPI tasks.
integer nwauxio
Number of auxiliary variables that are only included in output.
logical, dimension(:), allocatable w_write
if true write the w variable in output
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
logical, dimension(:), allocatable writelevel
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
character(len=std_len) base_filename
Base file name for simulation output, which will be followed by a number.
double precision, dimension(^nd, 2) writespshift
domain percentage cut off shifted from each boundary when converting data
integer, parameter unitconvert
integer, dimension(:,:), allocatable node
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:55
procedure(sub_check_params), pointer phys_te_images
Definition mod_physics.t:88
Module with all the methods that users can customize in AMRVAC.
procedure(aux_output), pointer usr_aux_output
procedure(special_convert), pointer usr_special_convert
Pointer to a tree_node.
Definition mod_forest.t:6