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