MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_slice.t
Go to the documentation of this file.
1 !> Writes D-1 slice, can do so in various formats, depending on slice_type
2 module mod_slice
4  use mod_comm_lib, only: mpistop
5  implicit none
6 
7  !> Maximum number of slices
8  integer, parameter :: nslicemax=1000
9 
10  !> Slice coordinates, see @ref slices.md
11  double precision :: slicecoord(nslicemax)
12 
13  !> the file number of slices
14  integer :: slicenext
15 
16  !> Number of slices to output
17  integer :: nslices
18 
19  !> The slice direction for each slice
20  integer :: slicedir(nslicemax)
21 
22  !> choose data type of slice: vtu, vtuCC, dat, or csv
23  character(len=std_len) :: slice_type
24 
25  !> tag for MPI message
26  integer, private :: itag
27 
28 contains
29 
30  subroutine write_slice
32  ! Writes a D-1 slice
33  ! by Oliver Porth
34  ! 22.Nov 2011
35  integer :: islice
36 
37  do islice=1,nslices
38  call put_slice(slicedir(islice),slicecoord(islice))
39  end do
40 
42  end subroutine write_slice
43 
44  subroutine put_slice(dir,xslice)
47  ! Writes a D-1 slice
48  ! For ONED simulations, output will be appended to one csv-file per slice
49  ! slices are sensitive to the saveprim switch and
50  ! can contain auxiliary io variables (nwauxio)
51  ! Thus csv or vtu(CC)-files with primitive variables are obtained.
52  ! when slice_type='dat', we save a D-1 dat file for potential restarts
53  ! by Oliver Porth
54  ! 22.Nov 2011
55  integer, intent(in) :: dir
56  double precision, intent(in) :: xslice
57  ! .. local ..
58  integer :: Njgrid, jgrid
59  integer, dimension(ndim-1) :: ixsubGlo, ixsubGhi
60  integer, dimension(ndim-1) :: ixsubMlo, ixsubMhi
61  integer :: size_subblock_io, nx^D, slice_fh, nwexpand
62  integer :: type_subblock_io, type_subblockC_io, type_subblock_x_io, type_subblockC_x_io
63  integer, dimension(ndim) :: sizes, subsizes, start
64  double precision,dimension(0:nw+nwauxio) :: normconv
65 
66  ! Preamble:
67  nx^d=ixmhi^d-ixmlo^d+1;
68  slice_fh=unitslice
69 
70  if (ndim==1) then
71  nwexpand = nwauxio
72  else
73  if(slice_type/='dat')then
74  nwexpand = nwauxio
75  else
76  nwexpand = 0
77  endif
78  end if
79 
80  ! Do a last consistency check:
81  select case(dir)
82  {case(^d)
83  if(xslice<xprobmin^d.or.xslice>xprobmax^d) &
84  call mpistop("slice out of bounds")
85  \}
86  end select
87 
88  ! Traverse the forest and fill nodes:
89  call select_slice(dir,xslice,.false.,slice_fh,normconv)
90 
91  ! Create the MPI-datatype and select indices:
92  {^ifthreed
93  select case(dir)
94  case (1)
95  sizes(1) = ixghi2; sizes(2) = ixghi3;
96  ixsubglo(1) = ixglo2; ixsubglo(2) = ixglo3;
97  ixsubghi(1) = ixghi2; ixsubghi(2) = ixghi3;
98  subsizes(1)=nx2;subsizes(2)=nx3;
99  start(1)=ixmlo2-1;start(2)=ixmlo3-1;
100  size_subblock_io=nx2*nx3*(nw+nwexpand)*size_double
101  case (2)
102  sizes(1) = ixghi1; sizes(2) = ixghi3;
103  ixsubglo(1) = ixglo1; ixsubglo(2) = ixglo3;
104  ixsubghi(1) = ixghi1; ixsubghi(2) = ixghi3;
105  subsizes(1)=nx1;subsizes(2)=nx3;
106  start(1)=ixmlo1-1;start(2)=ixmlo3-1;
107  size_subblock_io=nx1*nx3*(nw+nwexpand)*size_double
108  case (3)
109  ixsubglo(1) = ixglo1; ixsubglo(2) = ixglo2;
110  ixsubghi(1) = ixghi1; ixsubghi(2) = ixghi2;
111  sizes(1) = ixghi1; sizes(2) = ixghi2;
112  subsizes(1)=nx1;subsizes(2)=nx2;
113  start(1)=ixmlo1-1;start(2)=ixmlo2-1;
114  size_subblock_io=nx1*nx2*(nw+nwexpand)*size_double
115  case default
116  call mpistop("slice direction not clear in put_slice")
117  end select
118  }
119  {^iftwod
120  select case(dir)
121  case (1)
122  ixsubglo(1) = ixglo2; ixsubghi(1) = ixghi2;
123  sizes(1) = ixghi2
124  subsizes(1)=nx2
125  start(1)=ixmlo2-1
126  size_subblock_io=nx2*(nw+nwexpand)*size_double
127  case (2)
128  ixsubglo(1) = ixglo1; ixsubghi(1) = ixghi1;
129  sizes(1) = ixghi1
130  subsizes(1)=nx1
131  start(1)=ixmlo1-1
132  size_subblock_io=nx1*(nw+nwexpand)*size_double
133  case default
134  call mpistop("slice direction not clear in put_slice")
135  end select
136  }
137  {^ifoned
138  size_subblock_io=(nw+nwexpand)*size_double
139  }
140 
141  {^nooned
142  {^de&ixsubmlo(^de-1) = ixsubglo(^de-1)+nghostcells;}
143  {^de&ixsubmhi(^de-1) = ixsubghi(^de-1)-nghostcells;}
144  }
145 
146  sizes(ndim)=nw+nwexpand
147  subsizes(ndim)=nw+nwexpand
148  start(ndim)=0
149 
150  ! Types for center variables:
151  call mpi_type_create_subarray(ndim,sizes,subsizes,start, &
152  mpi_order_fortran,mpi_double_precision, &
153  type_subblock_io,ierrmpi)
154  call mpi_type_commit(type_subblock_io,ierrmpi)
155 
156  sizes(ndim)=^nd
157  subsizes(ndim)=^nd
158  start(ndim)=0
159  call mpi_type_create_subarray(ndim,sizes,subsizes,start, &
160  mpi_order_fortran,mpi_double_precision, &
161  type_subblock_x_io,ierrmpi)
162  call mpi_type_commit(type_subblock_x_io,ierrmpi)
163 
164 
165  ! Types for corner variables:
166  subsizes(1:ndim-1) = subsizes(1:ndim-1) + 1
167  start(1:ndim-1) = start(1:ndim-1) - 1
168  sizes(ndim)=nw+nwexpand
169  subsizes(ndim)=nw+nwexpand
170  start(ndim)=0
171 
172  call mpi_type_create_subarray(ndim,sizes,subsizes,start, &
173  mpi_order_fortran,mpi_double_precision, &
174  type_subblockc_io,ierrmpi)
175  call mpi_type_commit(type_subblockc_io,ierrmpi)
176 
177  sizes(ndim)=^nd
178  subsizes(ndim)=^nd
179  start(ndim)=0
180  call mpi_type_create_subarray(ndim,sizes,subsizes,start, &
181  mpi_order_fortran,mpi_double_precision, &
182  type_subblockc_x_io,ierrmpi)
183  call mpi_type_commit(type_subblockc_x_io,ierrmpi)
184 
185 
186  ! local number of sub-grids:
188 
189  ! Now output using various schemes:
190  if (ndim==1) then
191  call put_slice_zerod
192  else
193  select case(slice_type)
194  case ('csv')
195  call put_slice_csv
196  case ('dat')
197  call put_slice_dat
198  case ('vtu', 'vtuCC')
199  call put_slice_vtu
200  end select
201  end if
202 
203  do jgrid=1,njgrid
204  call dealloc_subnode(jgrid)
205  end do
206 
207  call mpi_type_free(type_subblock_io,ierrmpi)
208  call mpi_type_free(type_subblock_x_io,ierrmpi)
209  call mpi_type_free(type_subblockc_io,ierrmpi)
210  call mpi_type_free(type_subblockc_x_io,ierrmpi)
211 
212 
213  contains
214 
215  subroutine put_slice_vtu
216 
217  use mod_calculate_xw
218  character(len=1024) :: filename, xlabel
219  character(len=79) :: xxlabel
220  logical :: fileopen
221  character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
222  character(len=1024) :: outfilehead
223  integer :: status(MPI_STATUS_SIZE), ipe
224 
225  if (mype==0) then
226 
227  inquire(slice_fh,opened=fileopen)
228  if(.not.fileopen)then
229  ! generate filename:
230  write(xlabel,"(D9.2)")xslice
231  xxlabel=trim(xlabel)
232  if(xslice>=zero)then
233  write(xxlabel(1:1),"(a)") "+"
234  endif
235  write(filename,"(a,i1.1,a,i4.4,a)") trim(base_filename)//'_d',dir,'_x'//trim(xxlabel)//'_n',slicenext,'.vtu'
236  open(slice_fh,file=filename,status='unknown',form='formatted')
237  end if
238  ! get and write the header:
239  call getheadernames(wnamei,xandwnamei,outfilehead)
240  ! generate xml header
241  write(slice_fh,'(a)')'<?xml version="1.0"?>'
242  write(slice_fh,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
243  if(type_endian==1)then
244  write(slice_fh,'(a)')' version="0.1" byte_order="LittleEndian">'
245  else
246  write(slice_fh,'(a)')' version="0.1" byte_order="BigEndian">'
247  endif
248  write(slice_fh,'(a)')' <UnstructuredGrid>'
249  write(slice_fh,'(a)')'<FieldData>'
250  write(slice_fh,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
251  'NumberOfTuples="1" format="ascii">'
252  write(slice_fh,*) real(global_time*time_convert_factor)
253  write(slice_fh,'(a)')'</DataArray>'
254  write(slice_fh,'(a)')'</FieldData>'
255 
256  ! write to file:
257  do jgrid=1, njgrid
258  call write_slice_vtk(jgrid,slice_fh,wnamei)
259  end do
260 
261  ! create a recv buffer using allocate, will be deallocated at the end of the routine:
262  call alloc_subnode(njgrid+1,dir,nwauxio)
263 
264  end if
265 
266  ! Also communicate the normconv array since processor zero might not have it yet:
267  if (npe>1) then
268  do ipe=1,npe-1
269  do jgrid=1,morton_sub_stop(ipe)-morton_sub_start(ipe)+1
270  itag=morton_sub_start(ipe)+jgrid-1
271  itag=itag*5
272  if (ipe == mype ) then
273  call mpi_send(ps_sub(jgrid)%x,1,type_subblock_x_io,0,itag,icomm,ierrmpi)
274  call mpi_send(ps_sub(jgrid)%w,1,type_subblock_io,0,itag+1,icomm,ierrmpi)
275  call mpi_send(ps_sub(jgrid)%xC,1,type_subblockc_x_io,0,itag+2,icomm,ierrmpi)
276  call mpi_send(ps_sub(jgrid)%wC,1,type_subblockc_io,0,itag+3,icomm,ierrmpi)
277  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision,0,itag+4,icomm,ierrmpi)
278  end if
279  if (mype == 0) then
280  call mpi_recv(ps_sub(njgrid+1)%x,1,type_subblock_x_io,ipe,itag,icomm,status,ierrmpi)
281  call mpi_recv(ps_sub(njgrid+1)%w,1,type_subblock_io,ipe,itag+1,icomm,status,ierrmpi)
282  call mpi_recv(ps_sub(njgrid+1)%xC,1,type_subblockc_x_io,ipe,itag+2,icomm,status,ierrmpi)
283  call mpi_recv(ps_sub(njgrid+1)%wC,1,type_subblockc_io,ipe,itag+3,icomm,status,ierrmpi)
284  call mpi_recv(normconv,nw+nwauxio+1,mpi_double_precision,ipe,&
285  itag+4,icomm,status,ierrmpi)
286  call write_slice_vtk(njgrid+1,slice_fh,wnamei)
287  end if
288  end do
289  end do
290  endif
291 
292  if (mype==0) then
293 
294  write(slice_fh,'(a)')'</UnstructuredGrid>'
295  write(slice_fh,'(a)')'</VTKFile>'
296  close(slice_fh)
297  call dealloc_subnode(njgrid+1)
298 
299  end if
300 
301  end subroutine put_slice_vtu
302 
303  subroutine write_slice_vtk(jgrid,slice_fh,wnamei)
304 
305  ! this only works for 2D and 3D, 1D reduction (line to point) not allowed
306  integer, intent(in) :: jgrid, slice_fh
307  character(len=name_len), intent(in) :: wnamei(1:nw+nwauxio)
308  ! This remainder part only for more than 1D, but nesting with NOONED gives problems
309  {#IFNDEF D1
310  ! .. local ..
311  integer :: ixC^L, ixCC^L, nc, np, iw
312  integer :: nx^DM, nxC^DM, icell, ix^DM
313  double precision :: x_VTK(1:3)
314  integer :: VTK_type
315  double precision, parameter :: minvalue = 1.0d-99, maxvalue = 1.0d+99
316 
317  {^dm&ixccmin^dm = ixsubmlo(^dm);}
318  {^dm&ixccmax^dm = ixsubmhi(^dm);}
319  {^dm&ixcmin^dm = ixsubmlo(^dm)-1;}
320  {^dm&ixcmax^dm = ixsubmhi(^dm);}
321 
322  nx^dm=ixccmax^dm-ixccmin^dm+1;
323  nxc^dm=nx^dm+1;
324  nc={nx^dm*} ! Number of cells per subgrid
325  np={nxc^dm*} ! Number of corner points per subgrid
326 
327  ! we write out every grid as one VTK PIECE
328  write(slice_fh,'(a,i7,a,i7,a)') &
329  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
330 
331  !==============================
332  ! celldata or pointdata?
333  !==============================
334  select case(slice_type)
335 
336  case('vtuCC') ! celldata
337  write(slice_fh,'(a)')'<CellData>'
338  do iw=1,nw+nwauxio
339  if(iw<=nw) then
340  if(.not.w_write(iw)) cycle
341  endif
342  write(slice_fh,'(a,a,a)')&
343  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
344  write(slice_fh,'(200(1pe14.6))') {^dm&(|}roundoff_minmax(ps_sub(jgrid)%w(ix^dm,iw)*normconv(iw),minvalue,maxvalue),{ix^dm=ixccmin^dm,ixccmax^dm)}
345  write(slice_fh,'(a)')'</DataArray>'
346  enddo
347  write(slice_fh,'(a)')'</CellData>'
348 
349 
350  case('vtu') ! pointdata
351  write(slice_fh,'(a)')'<PointData>'
352  do iw=1,nw+nwauxio
353  if(iw<=nw) then
354  if(.not.w_write(iw)) cycle
355  endif
356  write(slice_fh,'(a,a,a)')&
357  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
358  write(slice_fh,'(200(1pe14.6))') {^dm&(|}roundoff_minmax(ps_sub(jgrid)%wC(ix^dm,iw)*normconv(iw),minvalue,maxvalue),{ix^dm=ixcmin^dm,ixcmax^dm)}
359  write(slice_fh,'(a)')'</DataArray>'
360  enddo
361  write(slice_fh,'(a)')'</PointData>'
362 
363 
364  end select
365  !==============================
366  ! Done: celldata or pointdata?
367  !==============================
368 
369  !==============================
370  ! output Cornerpoints
371  !==============================
372  write(slice_fh,'(a)')'<Points>'
373  write(slice_fh,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
374  ! write cell corner coordinates in a backward dimensional loop, always 3D output
375  {do ix^dmb=ixcmin^dmb,ixcmax^dmb \}
376  x_vtk(1:3)=zero;
377  x_vtk(1:ndim)=ps_sub(jgrid)%xC(ix^dm,1:ndim)*normconv(0);
378  write(slice_fh,'(3(1pe14.6))') x_vtk
379  {^dm&end do \}
380  write(slice_fh,'(a)')'</DataArray>'
381  write(slice_fh,'(a)')'</Points>'
382  !==============================
383  ! Done: output Cornerpoints
384  !==============================
385 
386  !==============================
387  ! cell Metainformation
388  !==============================
389  write(slice_fh,'(a)')'<Cells>'
390 
391  ! connectivity part
392  write(slice_fh,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
393 
394  {^dm& do ix^dmb=1,nx^dmb\}
395  {^iftwod
396  write(slice_fh,'(2(i7,1x))')ix1-1,ix1
397  }{^ifthreed
398  write(slice_fh,'(4(i7,1x))')(ix2-1)*nxc1+ix1-1, &
399  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
400  }{^dm& end do\}
401 
402  write(slice_fh,'(a)')'</DataArray>'
403 
404  ! offsets data array
405  write(slice_fh,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
406  do icell=1,nc
407  write(slice_fh,'(i7)') icell*(2**(^nd-1))
408  end do
409  write(slice_fh,'(a)')'</DataArray>'
410 
411  ! VTK cell type data array
412  write(slice_fh,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
413  ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
414  {^iftwod vtk_type=3 \}
415  {^ifthreed vtk_type=8 \}
416  do icell=1,nc
417  write(slice_fh,'(i2)') vtk_type
418  enddo
419  write(slice_fh,'(a)')'</DataArray>'
420 
421  write(slice_fh,'(a)')'</Cells>'
422  !==============================
423  ! Done: cell Metainformation
424  !==============================
425  write(slice_fh,'(a)')'</Piece>'
426 
427  }
428  end subroutine write_slice_vtk
429 
430  subroutine put_slice_csv
431 
432  use mod_calculate_xw
433  character(len=1024) :: filename, xlabel
434  character(len=79) :: xxlabel
435  logical :: fileopen
436  character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
437  character(len=1024) :: outfilehead
438  integer :: iw, ipe, itag
439  character(len=1024) :: line
440  integer :: status(MPI_STATUS_SIZE)
441 
442  if (mype==0) then
443  inquire(slice_fh,opened=fileopen)
444  if(.not.fileopen)then
445  ! generate filename:
446  write(xlabel,"(D9.2)")xslice
447  xxlabel=trim(xlabel)
448  if(xslice>=zero)then
449  write(xxlabel(1:1),"(a)") "+"
450  endif
451  write(filename,"(a,i1.1,a,i4.4,a)") trim(base_filename)//'_d',dir,'_x'//trim(xxlabel)//'_n',slicenext,'.csv'
452  open(slice_fh,file=filename,status='unknown',form='formatted')
453  end if
454  ! get and write the header:
455  call getheadernames(wnamei,xandwnamei,outfilehead)
456  line=''
457  do iw=1,ndim+nw+nwauxio-1
458  line = trim(line)//trim(xandwnamei(iw))//', '
459  end do
460  line = trim(line)//trim(xandwnamei(ndim+nw+nwauxio))
461  write(slice_fh,'(a)')trim(line)
462  ! create a recv buffer using allocate, will be deallocated at the end of the routine:
463  call alloc_subnode(njgrid+1,dir,nwauxio)
464 
465  ! write to file:
466  do jgrid=1, njgrid
467  call put_slice_line(jgrid,slice_fh)
468  end do
469  end if
470 
471  ! Also communicate the normconv array since processor zero might not have it yet:
472  if (npe>1) then
473  do ipe=1,npe-1
474  do jgrid=1,morton_sub_stop(ipe)-morton_sub_start(ipe)+1
475  itag=morton_sub_start(ipe)+jgrid-1
476  if (ipe == mype ) then
477  call mpi_send(ps_sub(jgrid)%x,1,type_subblock_x_io,0,itag,icomm,ierrmpi)
478  call mpi_send(ps_sub(jgrid)%w,1,type_subblock_io,0,itag,icomm,ierrmpi)
479  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision,0,itag,icomm,ierrmpi)
480  end if
481  if (mype == 0) then
482  call mpi_recv(ps_sub(njgrid+1)%x,1,type_subblock_x_io,ipe,itag,icomm,status,ierrmpi)
483  call mpi_recv(ps_sub(njgrid+1)%w,1,type_subblock_io,ipe,itag,icomm,status,ierrmpi)
484  call mpi_recv(normconv,nw+nwauxio+1,mpi_double_precision,ipe,&
485  itag,icomm,status,ierrmpi)
486  call put_slice_line(njgrid+1,slice_fh)
487  end if
488  end do
489  end do
490  endif
491 
492  if (mype==0) then
493  close(slice_fh)
494  call dealloc_subnode(njgrid+1)
495  end if
496 
497  end subroutine put_slice_csv
498 
499  subroutine put_slice_line(jout,file_handle)
500  integer, intent(in) :: jout, file_handle
501  ! .. local ..
502  character(len=1024) ::line, data
503  integer :: ix^D,idir,iw
504  double precision, parameter :: minvalue = 1.0d-99, maxvalue = 1.0d+99
505 
506  {^ifthreed
507  do ix2=ixsubmlo(2),ixsubmhi(2)
508  do ix1=ixsubmlo(1),ixsubmhi(1)
509  \}
510  {^iftwod
511  do ix1=ixsubmlo(1),ixsubmhi(1)
512  }
513  ! Format the line:
514  line = ''
515  do idir=1,ndim
516  {^ifthreed
517  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%x(ix1,ix2,idir),minvalue,maxvalue)
518  }
519  {^iftwod
520  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%x(ix1,idir),minvalue,maxvalue)
521  }
522  {^ifoned
523  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%x(idir),minvalue,maxvalue)
524  }
525 
526  line = trim(line)//trim(data)//', '
527  end do
528  do iw = 1,nw+nwauxio-1
529  {^ifthreed
530  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(ix1,ix2,iw)*normconv(iw),minvalue,maxvalue)
531  }
532  {^iftwod
533  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(ix1,iw)*normconv(iw),minvalue,maxvalue)
534  }
535  {^ifoned
536  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(iw)*normconv(iw),minvalue,maxvalue)
537  }
538  line = trim(line)//trim(data)//', '
539  end do
540  {^ifthreed
541  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(ix1,ix2,nw+nwauxio)*normconv(nw+nwauxio),minvalue,maxvalue)
542  }
543  {^iftwod
544  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(ix1,nw+nwauxio)*normconv(nw+nwauxio),minvalue,maxvalue)
545  }
546  line = trim(line)//trim(data)
547  write(file_handle,'(a)')trim(line)
548  {^iftwod
549  end do
550  }
551  {^ifthreed
552  end do
553  end do
554  \}
555 
556  end subroutine put_slice_line
557 
558  subroutine put_slice_dat
559 
560  integer, dimension(max_blocks) :: iorequest
561  integer, dimension(MPI_STATUS_SIZE,max_blocks) :: iostatus
562  integer(kind=MPI_OFFSET_KIND) :: offset
563  integer :: nsubleafs
564  character(len=1024) :: filename, xlabel
565  character(len=79) :: xxlabel
566  integer :: amode, status(MPI_STATUS_SIZE), iwrite
567 
568  nsubleafs=morton_sub_stop(npe-1)
569  ! generate filename
570  write(xlabel,"(D9.2)")xslice
571  xxlabel=trim(xlabel)
572  if(xslice>=zero)then
573  write(xxlabel(1:1),"(a)") "+"
574  endif
575  write(filename,"(a,i1.1,a,i4.4,a)") trim(base_filename)//'_d',dir,'_x'//trim(xxlabel)//'_n',slicenext,'.dat'
576 
577  if(mype==0) then
578  open(unit=slice_fh,file=filename,status='replace')
579  close(unit=slice_fh)
580  end if
581 
582  amode=ior(mpi_mode_create,mpi_mode_wronly)
583  call mpi_file_open(icomm,filename,amode,mpi_info_null,slice_fh,ierrmpi)
584  iorequest=mpi_request_null
585  iwrite=0
586 
587  do jgrid=1,njgrid
588  iwrite=iwrite+1
589  offset=int(size_subblock_io,kind=mpi_offset_kind) &
590  *int(morton_sub_start(mype)+jgrid-2,kind=mpi_offset_kind)
591  call mpi_file_iwrite_at(slice_fh,offset,ps_sub(jgrid)%w,1,type_subblock_io, &
592  iorequest(iwrite),ierrmpi)
593  end do
594 
595  if (iwrite>0) call mpi_waitall(iwrite,iorequest,iostatus,ierrmpi)
596  call mpi_barrier(icomm, ierrmpi)
597  call mpi_file_close(slice_fh,ierrmpi)
598 
599  if (mype==0) then
600  amode=ior(mpi_mode_append,mpi_mode_wronly)
601  call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null, &
602  slice_fh,ierrmpi)
603 
604  call select_slice(dir,xslice,.true.,slice_fh,normconv)
605 
606  {call mpi_file_write(slice_fh,subsizes(^de-1),1,mpi_integer,status,ierrmpi)\}
607 ! call MPI_FILE_WRITE(slice_fh,eqpar,neqpar+nspecialpar, &
608 ! MPI_DOUBLE_PRECISION,status,ierrmpi)
609  call mpi_file_write(slice_fh,nsubleafs,1,mpi_integer,status,ierrmpi)
610  call mpi_file_write(slice_fh,levmax_sub,1,mpi_integer,status,ierrmpi)
611  call mpi_file_write(slice_fh,ndim-1,1,mpi_integer,status,ierrmpi)
612  call mpi_file_write(slice_fh,ndir,1,mpi_integer,status,ierrmpi)
613  call mpi_file_write(slice_fh,nw,1,mpi_integer,status,ierrmpi)
614 ! call MPI_FILE_WRITE(slice_fh,neqpar+nspecialpar,1,MPI_INTEGER,status,ierrmpi)
615  call mpi_file_write(slice_fh,it,1,mpi_integer,status,ierrmpi)
616  call mpi_file_write(slice_fh,global_time,1,mpi_double_precision,status,ierrmpi)
617 
618  call mpi_file_close(slice_fh,ierrmpi)
619  end if
620 
621  end subroutine put_slice_dat
622 
623  subroutine put_slice_zerod
624 
625  use mod_calculate_xw
626  integer:: iw
627  character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
628  character(len=1024) :: outfilehead
629  logical, save :: opened=.false.
630  character(len=1024) ::line, data
631  character(len=1024) :: filename, xlabel
632  character(len=79) :: xxlabel
633  integer :: amode, iwrite, status(MPI_STATUS_SIZE)
634 
635  {^ifoned
636  ! generate filename:
637  write(xlabel,"(D9.2)")xslice
638  xxlabel=trim(xlabel)
639  if(xslice>=zero)then
640  write(xxlabel(1:1),"(a)") "+"
641  endif
642  write(filename,"(a,i1.1,a,i4.4,a)") trim(base_filename)//'_d',dir,'_x'//trim(xxlabel)//'.csv'
643 
644  ! Open for header:
645  if (.not. opened .and. mype==0) then
646  amode=ior(mpi_mode_create,mpi_mode_wronly)
647  amode=ior(amode,mpi_mode_append)
648  call mpi_file_open(mpi_comm_self,filename,amode, &
649  mpi_info_null,slice_fh,ierrmpi)
650  ! Create the header:
651  call getheadernames(wnamei,xandwnamei,outfilehead)
652  line = 'time'
653  do iw=1,nw+nwauxio+ndim
654  line = trim(line)//', '//trim(xandwnamei(iw))
655  enddo
656  ! Write header:
657  call mpi_file_write(slice_fh,line,len_trim(line), &
658  mpi_character,status,ierrmpi)
659  call mpi_file_write(slice_fh,achar(10),1,mpi_character,status,ierrmpi)
660  call mpi_file_close(slice_fh, ierrmpi)
661  opened=.true.
662  endif ! opened
663  ! Now we let the processor, which holds the data, write - therefore barrier.
664  call mpi_barrier(icomm,ierrmpi)
665 
666  ! There should be only one processor holding the point:
667  if (njgrid > 0) then
668  amode=ior(mpi_mode_create,mpi_mode_wronly)
669  amode=ior(amode,mpi_mode_append)
670  call mpi_file_open(mpi_comm_self,filename,amode, &
671  mpi_info_null,slice_fh,ierrmpi)
672 
673  ! Format the line:
674  write(data,"(es14.6)")global_time
675  line = trim(data)
676  write(data,"(es14.6)")ps_sub(1)%x
677  line = trim(line)//', '//trim(data)
678  do iw = 1,nw+nwauxio
679  write(data,"(es14.6)")ps_sub(1)%w(iw)*normconv(iw)
680  line = trim(line)//', '//trim(data)
681  end do
682  !
683  call mpi_file_write(slice_fh,trim(line),len_trim(line), &
684  mpi_character,status,ierrmpi)
685  call mpi_file_write(slice_fh,achar(10),1,mpi_character,status,ierrmpi)
686  call mpi_file_close(slice_fh, ierrmpi)
687  end if ! Njgrid > 0
688  }
689  end subroutine put_slice_zerod
690 
691  end subroutine put_slice
692 
693  subroutine select_slice(dir,xslice,writeonly,file_handle,normconv)
696  integer, intent(in) :: dir
697  double precision, intent(in) :: xslice
698  integer, intent(in) :: file_handle
699  logical, intent(in) :: writeonly
700  double precision,dimension(0:nw+nwauxio),intent(out) :: normconv
701  ! .. local ..
702  integer :: ig^D, jgrid, slice_fh, ipe, mylevmax
703  integer, dimension(nlevelshi) :: igslice
704 
705  jgrid = 0
706  mylevmax = 0
707 
708  ! Find the global slice index for every level:
709  call get_igslice(dir,xslice,igslice)
710 
711  ! Traverse forest to find grids indicating slice:
712  {^ifthreed
713  select case(dir)
714  case (1)
715  ig1 = igslice(1)
716  do ig3=1,ng3(1)
717  do ig2=1,ng2(1)
718  call traverse_slice(tree_root(ig^d))
719  end do
720  end do
721  case (2)
722  ig2 = igslice(1)
723  do ig3=1,ng3(1)
724  do ig1=1,ng1(1)
725  call traverse_slice(tree_root(ig^d))
726  end do
727  end do
728  case (3)
729  ig3 = igslice(1)
730  do ig2=1,ng2(1)
731  do ig1=1,ng1(1)
732  call traverse_slice(tree_root(ig^d))
733  end do
734  end do
735  case default
736  call mpistop("slice direction not clear in select_slice")
737  end select
738  }
739  {^iftwod
740  select case(dir)
741  case (1)
742  ig1 = igslice(1)
743  do ig2=1,ng2(1)
744  call traverse_slice(tree_root(ig^d))
745  end do
746  case (2)
747  ig2 = igslice(1)
748  do ig1=1,ng1(1)
749  call traverse_slice(tree_root(ig^d))
750  end do
751  case default
752  call mpistop("slice direction not clear in select_slice")
753  end select
754  }
755  {^ifoned
756  ig1 = igslice(1)
757  call traverse_slice(tree_root(ig^d))
758  }
759 
760  if (.not.writeonly) then
761  ! Synchronize the levmax_sub for output (only rank 0 needs it):
762  levmax_sub = mylevmax
763  call mpi_allreduce(mpi_in_place,levmax_sub,1,mpi_integer,mpi_max,icomm,ierrmpi)
764 
765  ! Communicate the subgrid indices according to new Morton sub-sfc:
766  morton_sub_start(:) = 1
767  do ipe=0,npe-1
768  call mpi_gather(jgrid,1,mpi_integer,morton_sub_stop,1,mpi_integer,ipe,icomm,ierrmpi)
769  end do
770 
771  do ipe = 0, npe-2
773  morton_sub_stop(ipe+1) = morton_sub_stop(ipe)+morton_sub_stop(ipe+1)
774  end do
775  end if
776 
777  contains
778 
779  recursive subroutine traverse_slice(tree)
780  implicit none
781  type(tree_node_ptr) :: tree
782  integer :: ic^d
783  integer, dimension(MPI_STATUS_SIZE) :: status
784 
785  if (writeonly) then
786  call mpi_file_write(file_handle,tree%node%leaf,1,mpi_logical, &
787  status,ierrmpi)
788  end if
789 
790  if (tree%node%leaf) then
791  if (tree%node%ipe == mype.and..not.writeonly) then
792  mylevmax = max(mylevmax,tree%node%level)
793  call fill_subnode(tree%node%igrid,tree%node%active,jgrid,dir,xslice,normconv)
794  end if
795  return
796  end if
797  ! We are out for leaves now, continue for branches
798 
799  ! Get the correct child:
800  select case (dir)
801  {case (^d)
802  ic^d = igslice(tree%node%level+1) - 2 * tree%node%ig^d + 2
803  \}
804  case default
805  call mpistop("slice direction not clear in traverse_slice")
806  end select
807 
808  ! Recursively descend into the correct child branch:
809  {^ifthreed
810  select case(dir)
811  case (1)
812  do ic3=1,2
813  do ic2=1,2
814  call traverse_slice(tree%node%child(ic^d))
815  end do
816  end do
817  case (2)
818  do ic3=1,2
819  do ic1=1,2
820  call traverse_slice(tree%node%child(ic^d))
821  end do
822  end do
823  case (3)
824  do ic2=1,2
825  do ic1=1,2
826  call traverse_slice(tree%node%child(ic^d))
827  end do
828  end do
829  case default
830  call mpistop("slice direction not clear in traverse_slice")
831  end select
832  }
833  {^iftwod
834  select case(dir)
835  case (1)
836  do ic2=1,2
837  call traverse_slice(tree%node%child(ic^d))
838  end do
839  case (2)
840  do ic1=1,2
841  call traverse_slice(tree%node%child(ic^d))
842  end do
843  case default
844  call mpistop("slice direction not clear in traverse_slice")
845  end select
846  }
847  {^ifoned
848  call traverse_slice(tree%node%child(ic^d))
849  }
850 
851  end subroutine traverse_slice
852 
853  end subroutine select_slice
854 
855  subroutine fill_subnode(igrid,active,jgrid,dir,xslice,normconv)
857  use mod_calculate_xw
858  integer, intent(in) :: igrid, dir
859  integer, intent(inout) :: jgrid
860  logical, intent(in) :: active
861  double precision, intent(in) :: xslice
862  double precision,dimension(0:nw+nwauxio),intent(out) :: normconv
863  ! .. local ..
864  double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP, xC
865  double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP, xCC
866  double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
867  double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
868  double precision, dimension(ixG^T) :: x_save
869  integer :: ixslice, nwexpand, ixC^L, ixCC^L
870  logical :: mask(ixG^T)
871 
872  mask=.false.
873  ! increase grid-count:
874  jgrid=jgrid+1
875  ! Allocate subdim solution array:
876  if (ndim==1) then
877  nwexpand = nwauxio
878  else
879  if(slice_type/='dat')then
880  nwexpand = nwauxio
881  else
882  nwexpand = 0
883  endif
884  end if
885  call alloc_subnode(jgrid,dir,nwexpand)
886  call fill_subnode_info(igrid,jgrid,dir)
887 
888  mask(ixm^t)=.true.
889 
890  ! Now hunt for the index closest to the slice:
891  {^ifoned
892  ixslice = minloc(dabs(xslice-ps(igrid)%x(:,dir)),1,mask(:))
893  }
894  {^iftwod
895  select case (dir)
896  case (1)
897  ixslice = minloc(dabs(xslice-ps(igrid)%x(:,ixmlo2,dir)),1,mask(:,ixmlo2))
898  case (2)
899  ixslice = minloc(dabs(xslice-ps(igrid)%x(ixmlo1,:,dir)),1,mask(ixmlo1,:))
900  case default
901  call mpistop("slice direction not clear in fill_subnode")
902  end select
903  }
904  {^ifthreed
905  select case (dir)
906  case (1)
907  ixslice = minloc(dabs(xslice-ps(igrid)%x(:,ixmlo2,ixmlo3,dir)),1,mask(:,ixmlo2,ixmlo3))
908  case (2)
909  ixslice = minloc(dabs(xslice-ps(igrid)%x(ixmlo1,:,ixmlo3,dir)),1,mask(ixmlo1,:,ixmlo3))
910  case (3)
911  ixslice = minloc(dabs(xslice-ps(igrid)%x(ixmlo1,ixmlo2,:,dir)),1,mask(ixmlo1,ixmlo2,:))
912  case default
913  call mpistop("slice direction not clear in fill_subnode")
914  end select
915  }
916 
917  call calc_x(igrid,xc,xcc)
918  ! Set the coordinate to be exactly on the slice:
919  xc(:^d&,dir) = xslice
920  xcc(:^d&,dir) = xslice
921  call calc_grid(unitslice,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
922  ixc^l,ixcc^l,.true.)
923  ! CC stands for CellCenter
924  ! C stands for Corner
925 
926  ! Fill the subdimensional solution and position array:
927  {^ifoned
928  ps_sub(jgrid)%w(1:nw+nwexpand) = wcc_tmp(ixslice,1:nw+nwexpand)
929  ps_sub(jgrid)%x(1:ndim) = xcc_tmp(ixslice,1:ndim)
930  ps_sub(jgrid)%wC(1:nw+nwexpand) = wc_tmp(ixslice,1:nw+nwexpand)
931  ps_sub(jgrid)%xC(1:ndim) = xc_tmp(ixslice,1:ndim)
932  }
933  {^iftwod
934  select case (dir)
935  case (1)
936  ps_sub(jgrid)%w(ixccmin2:ixccmax2,1:nw+nwexpand) &
937  = wcc_tmp(ixslice,ixccmin2:ixccmax2,1:nw+nwexpand)
938  ps_sub(jgrid)%x(ixccmin2:ixccmax2,1:ndim) &
939  = xcc_tmp(ixslice,ixccmin2:ixccmax2,1:ndim)
940  ps_sub(jgrid)%wC(ixcmin2:ixcmax2,1:nw+nwexpand) &
941  = wc_tmp(ixslice,ixcmin2:ixcmax2,1:nw+nwexpand)
942  ps_sub(jgrid)%xC(ixcmin2:ixcmax2,1:ndim) &
943  = xc_tmp(ixslice,ixcmin2:ixcmax2,1:ndim)
944  case (2)
945  ps_sub(jgrid)%w(ixccmin1:ixccmax1,1:nw+nwexpand) &
946  = wcc_tmp(ixccmin1:ixccmax1,ixslice,1:nw+nwexpand)
947  ps_sub(jgrid)%x(ixccmin1:ixccmax1,1:ndim) &
948  = xcc_tmp(ixccmin1:ixccmax1,ixslice,1:ndim)
949  ps_sub(jgrid)%wC(ixcmin1:ixcmax1,1:nw+nwexpand) &
950  = wc_tmp(ixcmin1:ixcmax1,ixslice,1:nw+nwexpand)
951  ps_sub(jgrid)%xC(ixcmin1:ixcmax1,1:ndim) &
952  = xc_tmp(ixcmin1:ixcmax1,ixslice,1:ndim)
953  case default
954  call mpistop("slice direction not clear in fill_subnode")
955  end select
956  }
957  {^ifthreed
958  select case (dir)
959  case (1)
960  ps_sub(jgrid)%w(ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:nw+nwexpand) = &
961  wcc_tmp(ixslice,ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:nw+nwexpand)
962  ps_sub(jgrid)%x(ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:ndim) = &
963  xcc_tmp(ixslice,ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:ndim)
964  ps_sub(jgrid)%wC(ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:nw+nwexpand) = &
965  wc_tmp(ixslice,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:nw+nwexpand)
966  ps_sub(jgrid)%xC(ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:ndim) = &
967  xc_tmp(ixslice,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:ndim)
968  case (2)
969  ps_sub(jgrid)%w(ixccmin1:ixccmax1,ixccmin3:ixccmax3,1:nw+nwexpand) = &
970  wcc_tmp(ixccmin1:ixccmax1,ixslice,ixccmin3:ixccmax3,1:nw+nwexpand)
971  ps_sub(jgrid)%x(ixccmin1:ixccmax1,ixccmin3:ixccmax3,1:ndim) = &
972  xcc_tmp(ixccmin1:ixccmax1,ixslice,ixccmin3:ixccmax3,1:ndim)
973  ps_sub(jgrid)%wC(ixcmin1:ixcmax1,ixcmin3:ixcmax3,1:nw+nwexpand) = &
974  wc_tmp(ixcmin1:ixcmax1,ixslice,ixcmin3:ixcmax3,1:nw+nwexpand)
975  ps_sub(jgrid)%xC(ixcmin1:ixcmax1,ixcmin3:ixcmax3,1:ndim) = &
976  xc_tmp(ixcmin1:ixcmax1,ixslice,ixcmin3:ixcmax3,1:ndim)
977  case (3)
978  ps_sub(jgrid)%w(ixccmin1:ixccmax1,ixccmin2:ixccmax2,1:nw+nwexpand) = &
979  wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixslice,1:nw+nwexpand)
980  ps_sub(jgrid)%x(ixccmin1:ixccmax1,ixccmin2:ixccmax2,1:ndim) = &
981  xcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixslice,1:ndim)
982  ps_sub(jgrid)%wC(ixcmin1:ixcmax1,ixcmin2:ixcmax2,1:nw+nwexpand) = &
983  wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixslice,1:nw+nwexpand)
984  ps_sub(jgrid)%xC(ixcmin1:ixcmax1,ixcmin2:ixcmax2,1:ndim) = &
985  xc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixslice,1:ndim)
986  case default
987  call mpistop("slice direction not clear in fill_subnode")
988  end select
989  }
990 
991  end subroutine fill_subnode
992 
993  subroutine alloc_subnode(jgrid,dir,nwexpand)
995  integer, intent(in) :: jgrid, dir, nwexpand
996 
997  ! take care, what comes out is not necessarily a right handed system!
998  {^ifoned
999  allocate(ps_sub(jgrid)%w(1:nw+nwexpand),ps_sub(jgrid)%x(1:ndim))
1000  allocate(ps_sub(jgrid)%wC(1:nw+nwexpand),ps_sub(jgrid)%xC(1:ndim))
1001  }
1002  {^iftwod
1003  select case (dir)
1004  case (1)
1005  allocate(ps_sub(jgrid)%w(ixglo2:ixghi2,1:nw+nwexpand),&
1006  ps_sub(jgrid)%x(ixglo2:ixghi2,1:ndim), &
1007  ps_sub(jgrid)%wC(ixglo2:ixghi2,1:nw+nwexpand),&
1008  ps_sub(jgrid)%xC(ixglo2:ixghi2,1:ndim))
1009  case (2)
1010  allocate(ps_sub(jgrid)%w(ixglo1:ixghi1,1:nw+nwexpand),&
1011  ps_sub(jgrid)%x(ixglo1:ixghi1,1:ndim), &
1012  ps_sub(jgrid)%wC(ixglo1:ixghi1,1:nw+nwexpand),&
1013  ps_sub(jgrid)%xC(ixglo1:ixghi1,1:ndim))
1014  case default
1015  call mpistop("slice direction not clear in alloc_subnode")
1016  end select
1017  }
1018  {^ifthreed
1019  select case (dir)
1020  case (1)
1021  allocate(ps_sub(jgrid)%w(ixglo2:ixghi2,ixglo3:ixghi3,1:nw+nwexpand),&
1022  ps_sub(jgrid)%x(ixglo2:ixghi2,ixglo3:ixghi3,1:ndim), &
1023  ps_sub(jgrid)%wC(ixglo2:ixghi2,ixglo3:ixghi3,1:nw+nwexpand),&
1024  ps_sub(jgrid)%xC(ixglo2:ixghi2,ixglo3:ixghi3,1:ndim))
1025  case (2)
1026  allocate(ps_sub(jgrid)%w(ixglo1:ixghi1,ixglo3:ixghi3,1:nw+nwexpand),&
1027  ps_sub(jgrid)%x(ixglo1:ixghi1,ixglo3:ixghi3,1:ndim), &
1028  ps_sub(jgrid)%wC(ixglo1:ixghi1,ixglo3:ixghi3,1:nw+nwexpand),&
1029  ps_sub(jgrid)%xC(ixglo1:ixghi1,ixglo3:ixghi3,1:ndim))
1030  case (3)
1031  allocate(ps_sub(jgrid)%w(ixglo1:ixghi1,ixglo2:ixghi2,1:nw+nwexpand),&
1032  ps_sub(jgrid)%x(ixglo1:ixghi1,ixglo2:ixghi2,1:ndim), &
1033  ps_sub(jgrid)%wC(ixglo1:ixghi1,ixglo2:ixghi2,1:nw+nwexpand),&
1034  ps_sub(jgrid)%xC(ixglo1:ixghi1,ixglo2:ixghi2,1:ndim))
1035  case default
1036  call mpistop("slice direction not clear in alloc_subnode")
1037  end select
1038  }
1039  end subroutine alloc_subnode
1040 
1041  subroutine dealloc_subnode(jgrid)
1043  integer, intent(in) :: jgrid
1044 
1045  if (jgrid==0) then
1046  call mpistop("trying to delete a non-existing grid in dealloc_subnode")
1047  end if
1048 
1049  deallocate(ps_sub(jgrid)%w,ps_sub(jgrid)%x,ps_sub(jgrid)%wC,ps_sub(jgrid)%xC)
1050 
1051  ! reset the global node info:
1052  node_sub(:,jgrid)=0
1053  rnode_sub(:,jgrid)=zero
1054 
1055  end subroutine dealloc_subnode
1056 
1057  subroutine fill_subnode_info(igrid,jgrid,dir)
1059  integer, intent(in) :: igrid,jgrid,dir
1060 
1061  node_sub(plevel_,jgrid)=node(plevel_,igrid)
1062  {^ifthreed
1063  select case(dir)
1064  case (1)
1065  node_sub(pig1_,jgrid)=node(pig2_,igrid)
1066  node_sub(pig2_,jgrid)=node(pig3_,igrid)
1067  rnode_sub(rpdx1_,jgrid)=rnode(rpdx2_,igrid)
1068  rnode_sub(rpdx2_,jgrid)=rnode(rpdx3_,igrid)
1069  rnode_sub(rpxmin1_,jgrid)=rnode(rpxmin2_,igrid)
1070  rnode_sub(rpxmin2_,jgrid)=rnode(rpxmin3_,igrid)
1071  case (2)
1072  node_sub(pig1_,jgrid)=node(pig1_,igrid)
1073  node_sub(pig2_,jgrid)=node(pig3_,igrid)
1074  rnode_sub(rpdx1_,jgrid)=rnode(rpdx1_,igrid)
1075  rnode_sub(rpdx2_,jgrid)=rnode(rpdx3_,igrid)
1076  rnode_sub(rpxmin1_,jgrid)=rnode(rpxmin1_,igrid)
1077  rnode_sub(rpxmin2_,jgrid)=rnode(rpxmin3_,igrid)
1078  case (3)
1079  node_sub(pig1_,jgrid)=node(pig1_,igrid)
1080  node_sub(pig2_,jgrid)=node(pig2_,igrid)
1081  rnode_sub(rpdx1_,jgrid)=rnode(rpdx1_,igrid)
1082  rnode_sub(rpdx2_,jgrid)=rnode(rpdx2_,igrid)
1083  rnode_sub(rpxmin1_,jgrid)=rnode(rpxmin1_,igrid)
1084  rnode_sub(rpxmin2_,jgrid)=rnode(rpxmin2_,igrid)
1085  case default
1086  call mpistop("slice direction not clear in fill_subnode_info")
1087  end select
1088  }
1089  {^iftwod
1090  select case(dir)
1091  case (1)
1092  node_sub(pig1_,jgrid)=node(pig2_,igrid)
1093  rnode_sub(rpdx1_,jgrid)=rnode(rpdx2_,igrid)
1094  rnode_sub(rpxmin1_,jgrid)=rnode(rpxmin2_,igrid)
1095  case (2)
1096  node_sub(pig1_,jgrid)=node(pig1_,igrid)
1097  rnode_sub(rpdx1_,jgrid)=rnode(rpdx1_,igrid)
1098  rnode_sub(rpxmin1_,jgrid)=rnode(rpxmin1_,igrid)
1099  case default
1100  call mpistop("slice direction not clear in fill_subnode_info")
1101  end select
1102  }
1103  {^ifoned
1104  }
1105 
1106  end subroutine fill_subnode_info
1107 
1108  subroutine get_igslice(dir,x,igslice)
1110  integer, intent(in) :: dir
1111  double precision, intent(in) :: x
1112  integer, dimension(nlevelshi), intent(out) :: igslice
1113  ! .. local ..
1114  integer :: level
1115  double precision :: distance
1116  !double precision :: xsgrid(ndim,nlevelshi),qs(ndim),xmgrid(ndim,3),xnew,xlgrid(ndim,3)
1117  !integer :: nbefore
1118 
1119  if (x.ne.x) &
1120  call mpistop("get_igslice: your slice position is NaN!")
1121 
1122  select case (dir)
1123  {case (^d)
1124  select case (stretch_type(^d))
1125  case (stretch_none) ! This dimension is unstretched
1126  if (x<=xprobmin^d) then
1127  igslice=1
1128  else if (x>=xprobmax^d) then
1129  do level=1,refine_max_level
1130  igslice(level)=ng^d(level)
1131  end do
1132  else
1133  distance=x-xprobmin^d
1134  do level = 1, refine_max_level
1135  igslice(level) = ceiling(distance/dg^d(level))
1136  end do
1137  end if
1138 
1139  case (stretch_uni) ! Uniform stretching
1140 
1141  if (x<=xprobmin^d) then
1142  igslice=1
1143  else if (x>=xprobmax^d) then
1144  do level=1,refine_max_level
1145  igslice(level)=ng^d(level)
1146  end do
1147  else
1148  distance=x-xprobmin^d
1149  do level=1,refine_max_level
1150  igslice(level)= ceiling(dlog(distance/dxfirst(level,^d)*(qstretch(level,^d)-1.d0)+1.d0)&
1151  /dlog(qstretch(level,^d))/dble(block_nx^d))
1152  end do
1153  end if
1154 
1155  case (stretch_symm) ! Symmetric stretching
1156  ! symmetric stretch about 0.5*(xprobmin+xprobmax)
1157  if (x<=xprobmin^d) then
1158  igslice=1
1159  else if (x>=xprobmax^d) then
1160  do level=1,refine_max_level
1161  igslice(level)=ng^d(level)
1162  end do
1163  else if(x<xprobmin^d+xstretch^d) then
1164  distance=xprobmin^d+xstretch^d-x
1165  ! stretch to left from xprobmin+xstretch
1166  do level=1,refine_max_level
1167  igslice(level) = nstretchedblocks(level,^d)/2-int(dlog(distance*(qstretch(level,^d)-one)/&
1168  dxfirst(level,^d)+one)/dlog(qstretch(level,^d))/dble(block_nx^d))
1169  end do
1170  else if(x>xprobmax^d-xstretch^d) then
1171  distance=x-xprobmax^d+xstretch^d
1172  ! stretch to right from xprobmax-xstretch
1173  do level=1,refine_max_level
1174  igslice(level) = ceiling(dlog(distance*(qstretch(level,^d)-one)/&
1175  dxfirst(level,^d)+one)/dlog(qstretch(level,^d))/dble(block_nx^d))+ng^d(level)-nstretchedblocks(level,^d)/2
1176  end do
1177  else
1178  ! possible non-stretched central part
1179  distance=x-xprobmin^d-xstretch^d
1180  do level=1,refine_max_level
1181  igslice(level)=nstretchedblocks(level,^d)/2+ceiling(distance/dg^d(level))
1182  end do
1183  end if
1184  !xsgrid(^D,1)=half*(xprobmax^D-xprobmin^D)&
1185  ! /(half*domain_nx^D-half*nstretchedblocks_baselevel(^D)*block_nx^D &
1186  ! +(one-qstretch_baselevel(^D)**(half*nstretchedblocks_baselevel(^D)*block_nx^D)) &
1187  ! /(one-qstretch_baselevel(^D)))
1188  !xmgrid(^D,1)=xsgrid(^D,1)*(domain_nx^D-nstretchedblocks_baselevel(^D)*block_nx^D)
1189  !xlgrid(^D,1)=xsgrid(^D,1)
1190 
1191  !if (x .ge. xmgrid(^D,1)*half) then
1192  ! xnew=x-xmgrid(^D,1)*half
1193  ! do level=1,refine_max_level
1194  ! qs(^D)=qstretch_baselevel(^D)**(one/2.d0**(level-1))
1195  ! if (level .gt. 1) xsgrid(^D,level)=xsgrid(^D,level-1)/(one+qs(^D))
1196  ! nbefore=(domain_nx^D/block_nx^D-nstretchedblocks_baselevel(^D)/2)*2**(level-1)
1197  ! igslice(level)=floor((dlog(one-xnew/xsgrid(^D,level)*(one-qs(^D)))/dlog(qs(^D))-1.e-16)/block_nx^D)+1+nbefore
1198  ! if (x>=xprobmax^D) igslice(level)=domain_nx^D*2**(level-1)
1199  ! end do
1200  !else if(x .ge. -xmgrid(^D,1)*half) then
1201  ! xnew=x+xmgrid(^D,1)*half
1202  ! do level=1,refine_max_level
1203  ! nbefore=nstretchedblocks_baselevel(^D)/2*2**(level-1)
1204  ! if (level .gt. 1) xlgrid(^D,level)=xlgrid(^D,level-1)/2.d0
1205  ! igslice(level)=floor(((xnew-1.e-16)/xlgrid(^D,level))/block_nx^D)+1+nbefore
1206  ! end do
1207  !else
1208  ! xnew=xmgrid(^D,1)*half-x
1209  ! do level=1,refine_max_level
1210  ! qs(^D)=qstretch_baselevel(^D)**(one/2.d0**(level-1))
1211  ! if (level .gt. 1) xsgrid(^D,level)=xsgrid(^D,level-1)/(one+qs(^D))
1212  ! igslice(level)=nstretchedblocks_baselevel(^D)/2*2**(level-1) &
1213  ! -floor((dlog(one-xnew/xsgrid(^D,level)*(one-qs(^D)))/dlog(qs(^D)))/block_nx^D)
1214  ! if (x<=xprobmin^D) igslice(level)=1
1215  ! end do
1216  !end if
1217 
1218  case default
1219  call mpistop("stretch type not supported by get_igslice")
1220  end select
1221  \}
1222  case default
1223  call mpistop("slice direction not clear in get_igslice")
1224  end select
1225 
1226  end subroutine get_igslice
1227 
1228  double precision function roundoff_minmax(val,minval,maxval)
1229  implicit none
1230  double precision,intent(in) :: val, minval, maxval
1231 
1232  roundoff_minmax = val
1233 
1234  if (abs(roundoff_minmax) .le. minval) then
1235  roundoff_minmax = 0.0d0
1236  end if
1237 
1238  if (roundoff_minmax .gt. maxval) then
1239  roundoff_minmax = maxval
1240  else if (roundoff_minmax .lt. -maxval) then
1241  roundoff_minmax = -maxval
1242  end if
1243 
1244  ! Replace NaN with maxval (e.g. Paraview chokes on ASCII NaN):
1246 
1247  end function roundoff_minmax
1248 
1249 end module mod_slice
1250 
recursive subroutine traverse_slice(tree)
Definition: mod_slice.t:780
subroutine put_slice_dat
Definition: mod_slice.t:559
subroutine put_slice_zerod
Definition: mod_slice.t:624
subroutine put_slice_csv
Definition: mod_slice.t:431
subroutine write_slice_vtk(jgrid, slice_fh, wnamei)
Definition: mod_slice.t:304
subroutine put_slice_line(jout, file_handle)
Definition: mod_slice.t:500
subroutine put_slice_vtu
Definition: mod_slice.t:216
Module with basic data types used in amrvac.
integer, parameter size_double
Size (in bytes) of double precision real.
Handles computations for coordinates and variables in output.
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 getheadernames(wnamei, xandwnamei, outfilehead)
get all variables names
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.
Definition: mod_comm_lib.t:208
Module with basic grid data structures.
Definition: mod_forest.t:2
integer, dimension(:), allocatable, save morton_sub_start
Definition: mod_forest.t:67
integer, dimension(:), allocatable, save morton_sub_stop
Definition: mod_forest.t:67
type(tree_node_ptr), dimension(:^d &), allocatable, save tree_root
Pointers to the coarse grid.
Definition: mod_forest.t:29
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, parameter unitslice
integer, parameter stretch_uni
Unidirectional stretching from a side.
double precision global_time
The global simulation time.
double precision xstretch
physical extent of stretched border in symmetric stretching
double precision, dimension(:), allocatable dg
extent of grid blocks in domain per dimension, in array over levels
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision time_convert_factor
Conversion factor for time unit.
integer, dimension(:,:), allocatable nstretchedblocks
(even) number of (symmetrically) stretched blocks per level and dimension
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 block_nx
number of cells for each dimension in grid block excluding ghostcells
integer, parameter plevel_
integer, dimension(:), allocatable, parameter d
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
double precision, dimension(:,:), allocatable qstretch
Stretching factors and first cell size for each AMR level and dimension.
integer npe
The number of MPI tasks.
integer nwauxio
Number of auxiliary variables that are only included in output.
integer, parameter stretch_none
No stretching.
logical, dimension(:), allocatable w_write
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer nghostcells
Number of ghost cells surrounding a grid.
integer, dimension(ndim) stretch_type
What kind of stretching is used per dimension.
double precision, dimension(:,:), allocatable rnode_sub
character(len=std_len) base_filename
Base file name for simulation output, which will be followed by a number.
integer refine_max_level
Maximal number of AMR levels.
integer, parameter stretch_symm
Symmetric stretching around the center.
double precision, dimension(:,:), allocatable dxfirst
integer, dimension(:,:), allocatable node
integer, dimension(:,:), allocatable node_sub
Writes D-1 slice, can do so in various formats, depending on slice_type.
Definition: mod_slice.t:2
subroutine fill_subnode_info(igrid, jgrid, dir)
Definition: mod_slice.t:1058
subroutine get_igslice(dir, x, igslice)
Definition: mod_slice.t:1109
subroutine fill_subnode(igrid, active, jgrid, dir, xslice, normconv)
Definition: mod_slice.t:856
integer slicenext
the file number of slices
Definition: mod_slice.t:14
integer, parameter nslicemax
Maximum number of slices.
Definition: mod_slice.t:8
subroutine alloc_subnode(jgrid, dir, nwexpand)
Definition: mod_slice.t:994
subroutine put_slice(dir, xslice)
Definition: mod_slice.t:45
subroutine write_slice
Definition: mod_slice.t:31
character(len=std_len) slice_type
choose data type of slice: vtu, vtuCC, dat, or csv
Definition: mod_slice.t:23
integer, dimension(nslicemax) slicedir
The slice direction for each slice.
Definition: mod_slice.t:20
double precision function roundoff_minmax(val, minval, maxval)
Definition: mod_slice.t:1229
integer nslices
Number of slices to output.
Definition: mod_slice.t:17
double precision, dimension(nslicemax) slicecoord
Slice coordinates, see Slice output.
Definition: mod_slice.t:11
subroutine select_slice(dir, xslice, writeonly, file_handle, normconv)
Definition: mod_slice.t:694
subroutine dealloc_subnode(jgrid)
Definition: mod_slice.t:1042
Pointer to a tree_node.
Definition: mod_forest.t:6