1 #include "cpp_macros.h"
11 integer,
parameter :: af_dat_file_version = 3
16 type(
box_t),
intent(in) :: box
17 integer,
intent(in) :: n_var
19 new_vars(DTIMES(0:box%n_cell+1), n_var)
24 integer,
intent(in) :: my_unit
28 public :: af_write_tree
29 public :: af_read_tree
30 public :: af_tree_copy_variable
31 public :: af_write_vtk
32 public :: af_write_numpy
34 public :: af_write_plane
36 public :: af_write_silo
37 public :: af_write_line
42 subroutine af_write_tree(tree, filename, write_other_data)
43 type(
af_t),
intent(in) :: tree
44 character(len=*),
intent(in) :: filename
46 integer :: my_unit, lvl, id, n
47 character(len=400) :: fname
49 fname = trim(filename) //
".dat"
51 open(newunit=my_unit, file=trim(fname), form=
'unformatted', &
52 access=
'stream', status=
'replace')
54 write(my_unit) af_dat_file_version
55 write(my_unit) tree%ready
56 write(my_unit) tree%box_limit
57 write(my_unit) tree%highest_lvl
58 write(my_unit) tree%highest_id
59 write(my_unit) tree%n_cell
60 write(my_unit) tree%n_var_cell
61 write(my_unit) tree%n_var_face
62 write(my_unit) tree%coord_t
63 write(my_unit) tree%coarse_grid_size
64 write(my_unit) tree%periodic
65 write(my_unit) tree%r_base
66 write(my_unit) tree%dr_base
68 write(my_unit) tree%cc_names
69 write(my_unit) tree%fc_names
70 write(my_unit) tree%cc_num_copies
71 write(my_unit) tree%cc_write_output
72 write(my_unit) tree%cc_write_binary
73 write(my_unit) tree%fc_write_binary
75 write(my_unit) tree%n_removed_ids
76 write(my_unit) tree%removed_ids(1:tree%n_removed_ids)
80 do lvl = 1, tree%highest_lvl
81 write(my_unit)
size(tree%lvls(lvl)%ids)
82 write(my_unit) tree%lvls(lvl)%ids
84 write(my_unit)
size(tree%lvls(lvl)%leaves)
85 write(my_unit) tree%lvls(lvl)%leaves
87 write(my_unit)
size(tree%lvls(lvl)%parents)
88 write(my_unit) tree%lvls(lvl)%parents
91 do id = 1, tree%highest_id
92 associate(box => tree%boxes(id))
93 write(my_unit) box%in_use
94 if (.not. box%in_use) cycle
96 write(my_unit) box%n_cell
97 write(my_unit) box%n_bc
98 write(my_unit) box%n_stencils
99 write(my_unit) box%lvl
100 write(my_unit) box%tag
101 write(my_unit) box%ix
102 write(my_unit) box%parent
103 write(my_unit) box%children
104 write(my_unit) box%neighbors
105 write(my_unit) box%neighbor_mat
106 write(my_unit) box%dr
107 write(my_unit) box%r_min
108 write(my_unit) box%coord_t
110 do n = 1, tree%n_var_cell
111 if (tree%cc_write_binary(n))
then
112 write(my_unit) box%cc(dtimes(:), n)
116 do n = 1, tree%n_var_face
117 if (tree%fc_write_binary(n))
then
118 write(my_unit) box%fc(dtimes(:), :, n)
123 if (box%n_bc > 0)
then
124 write(my_unit) box%bc_index_to_nb
125 write(my_unit) box%nb_to_bc_index
126 write(my_unit) box%bc_type
127 write(my_unit) box%bc_val
128 write(my_unit) box%bc_coords
132 if (box%n_stencils > 0)
then
133 do n = 1, box%n_stencils
134 associate(stncl => box%stencils(n))
135 write(my_unit) stncl%key
136 write(my_unit) stncl%shape
137 write(my_unit) stncl%stype
138 write(my_unit) stncl%cylindrical_gradient
140 if (
allocated(stncl%c))
then
141 write(my_unit)
size(stncl%c)
142 write(my_unit) stncl%c
147 if (
allocated(stncl%v))
then
148 write(my_unit)
size(stncl%v, 1)
149 write(my_unit) stncl%v
154 if (
allocated(stncl%f))
then
156 write(my_unit) stncl%f
161 if (
allocated(stncl%bc_correction))
then
163 write(my_unit) stncl%bc_correction
168 if (
allocated(stncl%sparse_ix))
then
169 write(my_unit)
size(stncl%sparse_ix, 2)
170 write(my_unit) stncl%sparse_ix
175 if (
allocated(stncl%sparse_v))
then
176 write(my_unit)
size(stncl%sparse_v, 1)
177 write(my_unit) stncl%sparse_v
187 write(my_unit)
present(write_other_data)
189 if (
present(write_other_data))
then
190 call write_other_data(my_unit)
194 print *,
"af_write_tree: written " // trim(fname)
195 end subroutine af_write_tree
198 subroutine af_read_tree(tree, filename, read_other_data)
200 type(af_t),
intent(inout) :: tree
201 character(len=*),
intent(in) :: filename
203 integer :: my_unit, lvl, n, id, version, k, m, nc
204 logical :: other_data_present
206 open(newunit=my_unit, file=trim(filename), form=
'unformatted', &
207 access=
'stream', status=
'old', action=
'read')
209 read(my_unit) version
210 if (version /= af_dat_file_version)
then
211 print *,
"File version read:", version
212 print *,
"File version required:", af_dat_file_version
213 error stop
"af_read_tree: incompatible file versions"
216 read(my_unit) tree%ready
217 read(my_unit) tree%box_limit
218 read(my_unit) tree%highest_lvl
219 read(my_unit) tree%highest_id
220 read(my_unit) tree%n_cell
221 read(my_unit) tree%n_var_cell
222 read(my_unit) tree%n_var_face
223 read(my_unit) tree%coord_t
224 read(my_unit) tree%coarse_grid_size
225 read(my_unit) tree%periodic
226 read(my_unit) tree%r_base
227 read(my_unit) tree%dr_base
232 if (.not.
allocated(tree%cc_auto_vars)) &
233 allocate(tree%cc_auto_vars(0))
234 if (.not.
allocated(tree%cc_func_vars)) &
235 allocate(tree%cc_func_vars(0))
237 read(my_unit) tree%cc_names
238 read(my_unit) tree%fc_names
239 read(my_unit) tree%cc_num_copies
240 read(my_unit) tree%cc_write_output
241 read(my_unit) tree%cc_write_binary
242 read(my_unit) tree%fc_write_binary
244 allocate(tree%removed_ids(tree%box_limit))
245 read(my_unit) tree%n_removed_ids
246 read(my_unit) tree%removed_ids(1:tree%n_removed_ids)
248 do lvl = 1, tree%highest_lvl
250 allocate(tree%lvls(lvl)%ids(n))
251 read(my_unit) tree%lvls(lvl)%ids
254 allocate(tree%lvls(lvl)%leaves(n))
255 read(my_unit) tree%lvls(lvl)%leaves
258 allocate(tree%lvls(lvl)%parents(n))
259 read(my_unit) tree%lvls(lvl)%parents
262 do lvl = tree%highest_lvl+1, af_max_lvl
263 allocate(tree%lvls(lvl)%ids(0))
264 allocate(tree%lvls(lvl)%leaves(0))
265 allocate(tree%lvls(lvl)%parents(0))
268 allocate(tree%boxes(tree%box_limit))
270 do id = 1, tree%highest_id
271 associate(box => tree%boxes(id))
272 read(my_unit) box%in_use
273 if (.not. box%in_use) cycle
275 read(my_unit) box%n_cell
276 read(my_unit) box%n_bc
277 read(my_unit) box%n_stencils
278 read(my_unit) box%lvl
279 read(my_unit) box%tag
281 read(my_unit) box%parent
282 read(my_unit) box%children
283 read(my_unit) box%neighbors
284 read(my_unit) box%neighbor_mat
286 read(my_unit) box%r_min
287 read(my_unit) box%coord_t
289 call af_init_box(tree, id)
291 do n = 1, tree%n_var_cell
292 if (tree%cc_write_binary(n))
then
293 read(my_unit) box%cc(dtimes(:), n)
297 do n = 1, tree%n_var_face
298 if (tree%fc_write_binary(n))
then
299 read(my_unit) box%fc(dtimes(:), :, n)
304 if (box%n_bc > 0)
then
305 read(my_unit) box%bc_index_to_nb
306 read(my_unit) box%nb_to_bc_index
307 read(my_unit) box%bc_type
308 read(my_unit) box%bc_val
309 read(my_unit) box%bc_coords
313 if (box%n_stencils > 0)
then
314 allocate(box%stencils(box%n_stencils))
316 do n = 1, box%n_stencils
317 associate(stncl => box%stencils(n))
318 read(my_unit) stncl%key
319 read(my_unit) stncl%shape
320 read(my_unit) stncl%stype
321 read(my_unit) stncl%cylindrical_gradient
326 read(my_unit) stncl%c
331 allocate(stncl%v(k, dtimes(nc)))
332 read(my_unit) stncl%v
337 allocate(stncl%f(dtimes(nc)))
338 read(my_unit) stncl%f
343 allocate(stncl%bc_correction(dtimes(nc)))
344 read(my_unit) stncl%bc_correction
349 allocate(stncl%sparse_ix(ndim, k))
350 read(my_unit) stncl%sparse_ix
354 if (k > 0 .and. m > 0)
then
355 allocate(stncl%sparse_v(m, k))
356 read(my_unit) stncl%sparse_v
364 read(my_unit) other_data_present
366 if (
present(read_other_data))
then
367 if (other_data_present)
then
368 call read_other_data(my_unit)
370 error stop
"af_read_tree: other data is not present"
375 end subroutine af_read_tree
377 subroutine af_tree_copy_variable(tree_from, ivs_from, tree_to, ivs_to)
379 type(af_t),
intent(in) :: tree_from
380 integer,
intent(in) :: ivs_from(:)
381 type(af_t),
intent(inout) :: tree_to
382 integer,
intent(in) :: ivs_to(:)
383 integer :: lvl, id, n, nc, i, j, k
388 do lvl = 1, tree_to%highest_lvl
390 do n = 1,
size(tree_to%lvls(lvl)%leaves)
391 id = tree_to%lvls(lvl)%leaves(n)
392 nc = tree_to%boxes(id)%n_cell
395 rr = af_r_cc(tree_to%boxes(id), [ijk])
396 tree_to%boxes(id)%cc(ijk, ivs_to) = &
397 af_interp1(tree_from, rr, [ivs_from], success)
398 if (.not. success) error stop
"af_tree_copy_variable: interpolation error"
405 end subroutine af_tree_copy_variable
408 subroutine af_write_line(tree, filename, ivs, r_min, r_max, n_points)
410 type(af_t),
intent(in) :: tree
411 character(len=*),
intent(in) :: filename
412 integer,
intent(in) :: ivs(:)
413 real(dp),
intent(in) :: r_min(ndim)
414 real(dp),
intent(in) :: r_max(ndim)
415 integer,
intent(in) :: n_points
417 integer,
parameter :: my_unit = 100
418 character(len=400) :: fname
420 real(dp) :: r(ndim), dr_vec(ndim)
421 real(dp),
allocatable :: line_data(:, :)
425 dr_vec = (r_max - r_min) / max(1, n_points-1)
427 allocate(line_data(n_cc+ndim, n_points))
431 r = r_min + (i-1) * dr_vec
432 line_data(1:ndim, i) = r
433 line_data(ndim+1:ndim+n_cc, i) = af_interp1(tree, r, ivs, success)
434 if (.not. success) error stop
"af_write_line: interpolation error"
438 fname = trim(filename) //
".txt"
441 open(my_unit, file=trim(fname), action=
"write")
443 write(my_unit,
'(A)', advance=
"no")
"# x"
445 write(my_unit,
'(A)', advance=
"no")
"# x y"
447 write(my_unit,
'(A)', advance=
"no")
"# x y z"
450 write(my_unit,
'(A)', advance=
"no")
" "//trim(tree%cc_names(ivs(i)))
452 write(my_unit,
'(A)')
""
456 write(my_unit, *) line_data(:, i)
460 end subroutine af_write_line
466 subroutine af_write_plane(tree, filename, ivs, r_min, r_max, n_pixels)
468 type(af_t),
intent(in) :: tree
469 character(len=*),
intent(in) :: filename
470 integer,
intent(in) :: ivs(:)
471 real(dp),
intent(in) :: r_min(ndim)
472 real(dp),
intent(in) :: r_max(ndim)
473 integer,
intent(in) :: n_pixels(2)
475 integer,
parameter :: my_unit = 100
476 character(len=100) :: fmt_string
477 character(len=400) :: fname
478 integer :: i, j, n_cc, dim_unused, n_points(3)
479 real(dp) :: r(ndim), dvec(3)
480 real(dp) :: v1(ndim), v2(ndim)
481 real(dp),
allocatable :: pixel_data(:, :, :)
486 dvec(1:2) = r_max(1:2) - r_min(1:2)
489 n_points = [n_pixels(1), n_pixels(2), 1]
490 v1 = [dvec(1), 0.0_dp] / (n_pixels(1) - 1)
491 v2 = [0.0_dp, dvec(2)] / (n_pixels(2) - 1)
494 dim_unused = minloc(abs(dvec), 1)
496 select case (dim_unused)
498 v1 = [0.0_dp, dvec(2), 0.0_dp] / (n_pixels(1) - 1)
499 v2 = [0.0_dp, 0.0_dp, dvec(3)] / (n_pixels(2) - 1)
500 n_points = [1, n_pixels(1), n_pixels(2)]
502 v1 = [dvec(1), 0.0_dp, 0.0_dp] / (n_pixels(1) - 1)
503 v2 = [0.0_dp, 0.0_dp, dvec(3)] / (n_pixels(2) - 1)
504 n_points = [n_pixels(1), 1, n_pixels(2)]
506 v1 = [dvec(1), 0.0_dp, 0.0_dp] / (n_pixels(1) - 1)
507 v2 = [0.0_dp, dvec(2), 0.0_dp] / (n_pixels(2) - 1)
508 n_points = [n_pixels(1), n_pixels(2), 1]
512 allocate(pixel_data(n_cc, n_pixels(1), n_pixels(2)))
515 do j = 1, n_pixels(2)
516 do i = 1, n_pixels(1)
517 r = r_min + (i-1) * v1 + (j-1) * v2
518 pixel_data(:, i, j) = af_interp1(tree, r, ivs, success)
519 if (.not. success) error stop
"af_write_plane: interpolation error"
525 write(fmt_string,
'(A,I0,A)')
'(', n_pixels(1),
'E20.8)'
528 fname = trim(filename) //
".vtk"
530 open(my_unit, file=trim(fname), action=
"write")
531 write(my_unit,
'(A)')
"# vtk DataFile Version 2.0"
532 write(my_unit,
'(A)') trim(filename)
533 write(my_unit,
'(A)')
"ASCII"
534 write(my_unit,
'(A)')
"DATASET STRUCTURED_POINTS"
535 write(my_unit,
'(A,3I10)')
"DIMENSIONS ", n_points
537 write(my_unit,
'(A,3E20.8)')
"ORIGIN ", [r_min(1), r_min(2), 0.0_dp]
538 write(my_unit,
'(A,3E20.8)')
"SPACING ", &
539 [v1(1) + v2(1), v1(2) + v2(2), 0.0_dp]
541 write(my_unit,
'(A,3E20.8)')
"ORIGIN ", r_min
542 write(my_unit,
'(A,3E20.8)')
"SPACING ", v1 + v2
544 write(my_unit,
'(A,2I0)')
"POINT_DATA ", product(n_points)
546 write(my_unit,
'(A)')
"SCALARS " // &
547 trim(tree%cc_names(ivs(i))) //
" double 1"
548 write(my_unit,
'(A)')
"LOOKUP_TABLE default"
549 write(my_unit, trim(fmt_string)) pixel_data(i, :, :)
552 end subroutine af_write_plane
557 subroutine af_write_vtk(tree, filename, n_cycle, time, ixs_cc, &
561 type(af_t),
intent(in) :: tree
562 character(len=*),
intent(in) :: filename
563 integer,
intent(in),
optional :: n_cycle
564 real(dp),
intent(in),
optional :: time
565 integer,
intent(in),
optional :: ixs_cc(:)
567 character(len=*),
intent(in),
optional :: add_names(:)
569 integer :: lvl, bc, bn, n, n_cells, n_nodes
570 integer :: ig, ijk, id, n_ix, c_ix, n_grids
571 integer :: cell_ix, node_ix, n_cycle_val
572 integer :: n_cc, n_add
573 integer,
parameter :: n_ch = af_num_children
574 integer :: nodes_per_box, cells_per_box
576 real(dp),
allocatable :: coords(:), cc_vars(:,:)
577 integer,
allocatable :: offsets(:), connects(:)
578 integer,
allocatable :: cell_types(:), icc_val(:)
580 character(len=400) :: fname
581 character(len=100),
allocatable :: var_names(:)
582 real(dp),
allocatable :: cc(dtimes(:), :)
587 if (.not. tree%ready) error stop
"Tree not ready"
588 time_val = 0.0_dp;
if (
present(time)) time_val = time
589 n_cycle_val = 0;
if (
present(n_cycle)) n_cycle_val = n_cycle
590 n_add = 0;
if (
present(add_names)) n_add =
size(add_names)
592 if (
present(add_names) .neqv.
present(add_vars)) &
593 stop
"af_write_vtk: both arguments (add_names, add_vars) needed"
595 if (
present(ixs_cc))
then
596 if (maxval(ixs_cc) > tree%n_var_cell .or. &
597 minval(ixs_cc) < 1) stop
"af_write_vtk: wrong indices given (ixs_cc)"
598 allocate(icc_val(
size(ixs_cc)))
601 call get_output_vars(tree, icc_val)
606 allocate(var_names(n_cc+n_add))
607 var_names(1:n_cc) = tree%cc_names(icc_val)
609 if (
present(add_names))
then
610 var_names(n_cc+1:n_cc+n_add) = add_names(:)
615 nodes_per_box = bn**ndim
616 cells_per_box = bc**ndim
618 allocate(cc(dtimes(0:bc+1), n_cc + n_add))
621 do lvl = 1, tree%highest_lvl
622 n_grids = n_grids +
size(tree%lvls(lvl)%leaves)
624 n_nodes = nodes_per_box * n_grids
625 n_cells = cells_per_box * n_grids
627 allocate(coords(ndim * n_nodes))
628 allocate(cc_vars(n_cells, n_cc+n_add))
629 allocate(offsets(cells_per_box * n_grids))
630 allocate(cell_types(cells_per_box * n_grids))
631 allocate(connects(n_ch * cells_per_box * n_grids))
643 do lvl = 1, tree%highest_lvl
644 do n = 1,
size(tree%lvls(lvl)%leaves)
645 id = tree%lvls(lvl)%leaves(n)
647 cell_ix = (ig-1) * cells_per_box
648 node_ix = (ig-1) * nodes_per_box
651 cc(:, 1:n_cc) = tree%boxes(id)%cc(:, icc_val)
653 if (
present(add_vars))
then
654 call add_vars(tree%boxes(id), &
655 cc(:, n_cc+1:n_cc+n_add), n_add)
660 coords(n_ix) = tree%boxes(id)%r_min(1) + &
661 (i-1) * tree%boxes(id)%dr(1)
666 n_ix = node_ix + i - 1
668 cc_vars(c_ix, :) = cc(i, :)
669 offsets(c_ix) = af_num_children * c_ix
670 connects(n_ch*(c_ix-1)+1:n_ch*c_ix) = [n_ix, n_ix+1]
673 cc(:, :, 1:n_cc) = tree%boxes(id)%cc(:, :, icc_val)
675 if (
present(add_vars))
then
676 call add_vars(tree%boxes(id), &
677 cc(:, :, n_cc+1:n_cc+n_add), n_add)
682 n_ix = 2 * (node_ix + (j-1) * bn + i)
683 coords(n_ix-1:n_ix) = tree%boxes(id)%r_min + &
684 [i-1,j-1] * tree%boxes(id)%dr
691 n_ix = node_ix + (j-1) * bn + i - 1
692 c_ix = cell_ix + (j-1) * bc + i
693 cc_vars(c_ix, :) = cc(i, j, :)
694 offsets(c_ix) = af_num_children * c_ix
695 connects(n_ch*(c_ix-1)+1:n_ch*c_ix) = [n_ix, n_ix+1, n_ix+bn, n_ix+bn+1]
699 cc(:, :, :, 1:n_cc) = tree%boxes(id)%cc(:, :, :, icc_val)
701 if (
present(add_vars))
then
702 call add_vars(tree%boxes(id), &
703 cc(:, :, :, n_cc+1:n_cc+n_add), n_add)
709 n_ix = 3 * (node_ix + (k-1) * bn2 + (j-1) * bn + i)
710 coords(n_ix-2:n_ix) = tree%boxes(id)%r_min + &
711 [i-1,j-1,k-1] * tree%boxes(id)%dr
720 n_ix = node_ix + (k-1) * bn2 + &
722 c_ix = cell_ix + (k-1) * bc**2 + &
724 cc_vars(c_ix, :) = cc(i, j, k, :)
725 offsets(c_ix) = 8 * c_ix
726 connects(n_ch*(c_ix-1)+1:n_ch*c_ix) = &
727 [n_ix, n_ix+1, n_ix+bn, n_ix+bn+1, &
728 n_ix+bn2, n_ix+bn2+1, n_ix+bn2+bn, n_ix+bn2+bn+1]
736 fname = trim(filename) //
".vtu"
738 call vtk_ini_xml(vtkf, trim(fname),
'UnstructuredGrid')
739 call vtk_dat_xml(vtkf,
"UnstructuredGrid", .true.)
740 call vtk_unstr_geo_xml(vtkf, coords, n_nodes, n_cells, ndim, n_cycle_val, time_val)
741 call vtk_unstr_con_xml(vtkf, connects, offsets, cell_types, n_cells)
742 call vtk_dat_xml(vtkf,
"CellData", .true.)
744 do n = 1, n_cc + n_add
745 call vtk_var_r8_xml(vtkf, trim(var_names(n)), cc_vars(:, n), n_cells)
748 call vtk_dat_xml(vtkf,
"CellData", .false.)
749 call vtk_unstr_geo_xml_close(vtkf)
750 call vtk_dat_xml(vtkf,
"UnstructuredGrid", .false.)
751 call vtk_end_xml(vtkf)
752 print *,
"af_write_vtk: written " // trim(fname)
753 end subroutine af_write_vtk
757 subroutine af_write_numpy(tree, filename, r_min, r_max, n_points, &
758 n_cycle, time, ixs_cc)
762 type(af_t),
intent(inout) :: tree
764 character(len=*),
intent(in) :: filename
765 integer,
intent(in),
optional :: n_points(ndim)
766 real(dp),
intent(in),
optional :: r_min(ndim)
767 real(dp),
intent(in),
optional :: r_max(ndim)
768 integer,
intent(in),
optional :: n_cycle
769 real(dp),
intent(in),
optional :: time
770 integer,
intent(in),
optional :: ixs_cc(:)
772 integer :: n, ijk, id, lvl, id_guess, nx(ndim)
773 integer :: n_cycle_val, n_cc, nc
775 real(dp) :: time_val, dr(ndim), r(ndim)
776 real(dp) :: rmin(ndim), rmax(ndim)
777 integer,
allocatable :: icc_val(:)
778 character(len=af_nlen),
allocatable :: var_names(:)
779 character(len=400) :: fname
780 real(dp),
allocatable :: cc(dtimes(:), :)
782 if (.not. tree%ready) error stop
"Tree not ready"
784 time_val = 0.0_dp;
if (
present(time)) time_val = time
785 n_cycle_val = 0;
if (
present(n_cycle)) n_cycle_val = n_cycle
787 if (
present(r_min)) rmin = r_min
788 rmax = tree%r_base + tree%dr_base * tree%coarse_grid_size
789 if (
present(r_max)) rmax = r_max
791 if (
present(ixs_cc))
then
792 if (maxval(ixs_cc) > tree%n_var_cell .or. &
793 minval(ixs_cc) < 1) stop
"af_write_vtk: wrong indices given (ixs_cc)"
794 allocate(icc_val(
size(ixs_cc)))
797 call get_output_vars(tree, icc_val)
801 if (
present(n_points))
then
805 outer:
do lvl = 1, tree%highest_lvl-1
806 do i = 1,
size(tree%lvls(lvl)%leaves)
807 id = tree%lvls(lvl)%leaves(i)
810 if (all(tree%boxes(id)%r_min <= rmax .and. &
811 tree%boxes(id)%r_min + nc * &
812 tree%boxes(id)%dr >= rmin))
exit outer
817 nx = nint((rmax - rmin)/af_lvl_dr(tree, lvl))
820 dr = (rmax - rmin)/nx
822 allocate(var_names(n_cc))
823 var_names(1:n_cc) = tree%cc_names(icc_val)
824 allocate(cc(dindex(nx), n_cc))
830 r = rmin + ([ijk] - 0.5_dp) * dr
831 cc(ijk, :) = af_interp1(tree, r, icc_val, success, id_guess)
832 if (.not. success) error stop
"af_write_plane: interpolation error"
837 r = rmin + ([ijk] - 0.5) * dr
838 cc(ijk, :) = af_interp1(tree, r, icc_val, success, id_guess)
839 if (.not. success) error stop
"af_write_plane: interpolation error"
846 r = rmin + ([ijk] - 0.5) * dr
847 cc(ijk, :) = af_interp1(tree, r, icc_val, success, id_guess)
848 if (.not. success) error stop
"af_write_plane: interpolation error"
855 n = len_trim(filename)
858 select case (filename(n-3:n))
862 call remove_file(filename)
866 fname = filename(:n-4)//trim(var_names(n))//
'.npy'
867 call save_npy(fname, cc(dtimes(:), n))
868 call add_to_zip(filename, fname, .false., var_names(n))
871 fname = filename(:n-4)//
'nx.npy'
873 call add_to_zip(filename, fname, .false.,
'nx')
875 fname = filename(:n-4)//
'r_min.npy'
877 call add_to_zip(filename, fname, .false.,
'r_min')
879 fname = filename(:n-4)//
'r_max.npy'
881 call add_to_zip(filename, fname, .false.,
'r_max')
883 fname = filename(:n-4)//
'dr.npy'
885 call add_to_zip(filename, fname, .false.,
'dr')
887 fname = filename(:n-4)//
'coord_t.npy'
888 call save_npy(fname, [tree%coord_t])
889 call add_to_zip(filename, fname, .false.,
'coord_t')
891 fname = filename(:n-4)//
'time_cycle.npy'
892 call save_npy(fname, [time_val, real(n_cycle_val, dp)])
893 call add_to_zip(filename, fname, .false.,
'time_cycle')
895 error stop
"Unknown file extension"
898 print *,
"af_write_numpy: written " // trim(filename)
899 end subroutine af_write_numpy
904 subroutine af_write_silo(tree, filename, n_cycle, time, ixs_cc, &
905 add_vars, add_names, max_lvl)
909 type(af_t),
intent(in) :: tree
910 character(len=*) :: filename
911 integer,
intent(in),
optional :: n_cycle
912 real(dp),
intent(in),
optional :: time
913 integer,
intent(in),
optional :: ixs_cc(:)
915 character(len=*),
intent(in),
optional :: add_names(:)
917 integer,
intent(in),
optional :: max_lvl
918 integer :: n, i, id, ix, lvl, n_boxes
919 integer :: dbix, n_cc, n_leaves, nc
920 integer :: highest_lvl
921 integer,
allocatable :: icc_val(:), id_leaves(:)
922 integer,
allocatable :: id_sorted_ix(:)
923 real(dp),
allocatable :: xmin_leaves(:)
924 real(dp),
allocatable :: xdata(:), ydata(:, :)
925 character(len=af_nlen) :: yname
926 character(len=400) :: fname
928 highest_lvl = tree%highest_lvl
929 if (
present(max_lvl)) highest_lvl = min(max_lvl, tree%highest_lvl)
931 if (
present(ixs_cc))
then
932 if (maxval(ixs_cc) > tree%n_var_cell .or. &
933 minval(ixs_cc) < 1) stop
"af_write_silo: wrong indices given (ixs_cc)"
934 allocate(icc_val(
size(ixs_cc)))
937 call get_output_vars(tree, icc_val)
940 if (
present(add_vars) .or.
present(add_names)) &
941 error stop
"add_vars / add_names not implemented in 1D"
948 do lvl = 1, highest_lvl
949 if (lvl < highest_lvl)
then
950 n = n +
size(tree%lvls(lvl)%leaves)
952 n = n +
size(tree%lvls(lvl)%ids)
957 allocate(ydata(n_leaves * nc, n_cc))
958 allocate(xdata(n_leaves * nc))
959 allocate(id_leaves(n_leaves))
960 allocate(xmin_leaves(n_leaves))
961 allocate(id_sorted_ix(n_leaves))
965 do lvl = 1, highest_lvl
966 if (lvl < highest_lvl)
then
967 n_boxes =
size(tree%lvls(lvl)%leaves)
968 id_leaves(n+1:n+n_boxes) = tree%lvls(lvl)%leaves
970 n_boxes =
size(tree%lvls(lvl)%ids)
971 id_leaves(n+1:n+n_boxes) = tree%lvls(lvl)%ids
977 xmin_leaves = tree%boxes(id_leaves)%r_min(1)
980 call mrgrnk(xmin_leaves, id_sorted_ix)
984 id = id_leaves(id_sorted_ix(n))
985 ydata(ix+1:ix+nc, :) = tree%boxes(id)%cc(1:nc, icc_val)
987 xdata(ix+i:ix+i) = af_r_cc(tree%boxes(id), [i])
991 fname = trim(filename) //
".silo"
992 call silo_create_file(trim(fname), dbix)
993 call silo_set_time_varying(dbix)
996 yname = tree%cc_names(icc_val(n))
997 call silo_add_curve(dbix, yname, xdata, ydata(:, n))
1000 call silo_close_file(dbix)
1001 print *,
"af_write_silo: written " // trim(fname)
1002 end subroutine af_write_silo
1006 subroutine af_write_silo(tree, filename, n_cycle, time, ixs_cc, &
1007 add_vars, add_names, add_curve_names, add_curve_dat, max_lvl)
1010 type(af_t),
intent(in) :: tree
1011 character(len=*) :: filename
1012 integer,
intent(in),
optional :: n_cycle
1013 real(dp),
intent(in),
optional :: time
1014 integer,
intent(in),
optional :: ixs_cc(:)
1016 character(len=*),
intent(in),
optional :: add_names(:), add_curve_names(:)
1017 real(dp),
intent(in),
optional :: add_curve_dat(:, :, :)
1020 integer,
intent(in),
optional :: max_lvl
1022 character(len=*),
parameter :: grid_name =
"gg", block_prefix =
"blk_"
1023 character(len=*),
parameter :: amr_name =
"mesh", meshdir =
"data"
1024 character(len=100),
allocatable :: grid_list(:), grid_list_block(:)
1025 character(len=100),
allocatable :: var_list(:, :), var_names(:)
1026 character(len=400) :: fname
1027 integer :: lvl, i, id, i_grid, iv, nc, n_grids_max
1028 integer :: n_cc, n_add, dbix, highest_lvl
1029 integer :: nx, nx_prev, ix
1030 integer :: n_cycle_val, ncurve, ncurve_add
1031 integer :: lo(ndim), hi(ndim), vlo(ndim), vhi(ndim)
1032 integer :: blo(ndim), bhi(ndim)
1033 logical :: lo_bnd(ndim), hi_bnd(ndim)
1034 integer,
allocatable :: ids(:), nb_ids(:), icc_val(:)
1035 logical,
allocatable :: box_done(:)
1036 real(dp) :: dr(ndim), r_min(ndim), time_val
1037 integer,
allocatable :: box_list(dtimes(:)), new_box_list(dtimes(:))
1038 real(dp),
allocatable :: var_data(dtimes(:),:), cc(dtimes(:), :)
1040 integer :: ny, ny_prev, iy
1043 integer :: nz, nz_prev, iz
1046 if (.not. tree%ready) stop
"Tree not ready"
1047 time_val = 0.0_dp;
if (
present(time)) time_val = time
1048 n_cycle_val = 0;
if (
present(n_cycle)) n_cycle_val = n_cycle
1049 n_add = 0;
if (
present(add_names)) n_add =
size(add_names)
1051 if (
present(add_names) .neqv.
present(add_vars)) &
1052 stop
"af_write_silo: both arguments (add_names, add_vars) needed"
1054 if (
present(add_curve_names) .neqv.
present(add_curve_dat)) &
1055 stop
"af_write_silo: both arguments (add_curve_names, add_curve_dat) needed"
1057 if (
present(add_curve_names))
then
1058 if (
size(add_curve_names, 1) .ne.
size(add_curve_dat, 1)) &
1059 stop
"af_write_silo: number of curve names and data do not agree"
1062 if (
present(ixs_cc))
then
1063 if (maxval(ixs_cc) > tree%n_var_cell .or. &
1064 minval(ixs_cc) < 1) stop
"af_write_silo: wrong indices given (ixs_cc)"
1065 allocate(icc_val(
size(ixs_cc)))
1068 call get_output_vars(tree, icc_val)
1071 highest_lvl = tree%highest_lvl
1072 if (
present(max_lvl)) highest_lvl = min(max_lvl, tree%highest_lvl)
1074 n_cc =
size(icc_val)
1076 allocate(var_names(n_cc+n_add))
1077 var_names(1:n_cc) = tree%cc_names(icc_val)
1079 if (
present(add_names))
then
1080 var_names(n_cc+1:n_cc+n_add) = add_names(:)
1085 do lvl = 1, highest_lvl
1086 if (lvl < highest_lvl)
then
1087 n_grids_max = n_grids_max +
size(tree%lvls(lvl)%leaves)
1089 n_grids_max = n_grids_max +
size(tree%lvls(lvl)%ids)
1093 allocate(grid_list(n_grids_max))
1094 allocate(grid_list_block(n_grids_max))
1095 allocate(var_list(n_cc+n_add, n_grids_max))
1096 allocate(box_done(tree%highest_id))
1099 allocate(cc(dtimes(0:nc+1), n_cc + n_add))
1101 fname = trim(filename) //
".silo"
1102 call silo_create_file(trim(fname), dbix)
1103 call silo_set_time_varying(dbix)
1104 call silo_mkdir(dbix, meshdir)
1107 if (
present(add_curve_names))
then
1108 ncurve_add =
size(add_curve_dat, 1)
1109 do ncurve = 1, ncurve_add
1110 call silo_add_curve(dbix, add_curve_names(ncurve), &
1111 add_curve_dat(ncurve, 1, :), add_curve_dat(ncurve, 2, :))
1117 do lvl = 1, highest_lvl
1118 do i = 1,
size(tree%lvls(lvl)%ids)
1119 id = tree%lvls(lvl)%ids(i)
1121 if (box_done(id)) cycle
1124 if (lvl < highest_lvl .and. af_has_children(tree%boxes(id))) cycle
1131 allocate(box_list(1))
1133 box_done(id) = .true.
1141 nb_ids = tree%boxes(ids)%neighbors(af_neighb_lowx)
1142 if (all(nb_ids > af_no_box))
then
1143 if (.not. any(box_done(nb_ids)) .and. &
1144 tree%boxes(nb_ids(1))%ix(1) < tree%boxes(ids(1))%ix(1) .and. &
1145 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1147 allocate(new_box_list(nx))
1148 new_box_list(1:1) = nb_ids
1149 new_box_list(2:) = box_list
1150 box_list = new_box_list
1151 box_done(nb_ids) = .true.
1152 deallocate(new_box_list)
1157 ids = box_list(nx:nx)
1158 nb_ids = tree%boxes(ids)%neighbors(af_neighb_highx)
1159 if (all(nb_ids > af_no_box))
then
1160 if (.not. any(box_done(nb_ids)) .and. &
1161 tree%boxes(nb_ids(1))%ix(1) > tree%boxes(ids(1))%ix(1) .and. &
1162 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1164 allocate(new_box_list(nx))
1165 new_box_list(nx:nx) = nb_ids
1166 new_box_list(1:nx-1) = box_list
1167 box_list = new_box_list
1168 box_done(nb_ids) = .true.
1169 deallocate(new_box_list)
1173 if (nx == nx_prev)
exit
1179 lo_bnd = af_is_phys_boundary(tree%boxes, id, af_low_neighbs)
1182 hi_bnd = af_is_phys_boundary(tree%boxes, id, af_high_neighbs)
1185 where (.not. lo_bnd) lo = lo - 1
1188 where (.not. hi_bnd) hi = hi + 1
1191 allocate(var_data(lo(1):hi(1), n_cc+n_add))
1196 cc(:, 1:n_cc) = tree%boxes(id)%cc(:, icc_val)
1197 if (
present(add_vars))
then
1198 call add_vars(tree%boxes(id), &
1199 cc(:, n_cc+1:n_cc+n_add), n_add)
1204 where ([ix] == 1 .and. .not. lo_bnd) blo = 0
1207 where ([ix] == [nx] .and. .not. hi_bnd) bhi = nc+1
1209 vlo = blo + ([ix]-1) * nc
1210 vhi = bhi + ([ix]-1) * nc
1212 var_data(vlo(1):vhi(1), :) = cc(blo(1):bhi(1), :)
1216 dr = tree%boxes(id)%dr
1217 r_min = tree%boxes(id)%r_min - (1 - lo) * dr
1219 write(grid_list(i_grid),
"(A,I0)") meshdir //
'/' // grid_name, i_grid
1220 call silo_add_grid(dbix, grid_list(i_grid), 1, &
1221 hi - lo + 2, r_min, dr, 1-lo, hi - [nx] * nc)
1222 write(grid_list_block(i_grid),
"(A,I0)") meshdir //
'/' // block_prefix &
1223 // grid_name, i_grid
1224 call silo_add_grid(dbix, grid_list_block(i_grid), 1, [nx+1], &
1225 tree%boxes(id)%r_min, nc*dr, [0], [0])
1227 do iv = 1, n_cc+n_add
1228 write(var_list(iv, i_grid),
"(A,I0)") meshdir //
'/' // &
1229 trim(var_names(iv)) //
"_", i_grid
1230 call silo_add_var(dbix, var_list(iv, i_grid), grid_list(i_grid), &
1231 pack(var_data(:, iv), .true.), hi-lo+1)
1234 deallocate(var_data)
1235 deallocate(box_list)
1237 allocate(box_list(1,1))
1239 box_done(id) = .true.
1248 ids = box_list(1, :)
1249 nb_ids = tree%boxes(ids)%neighbors(af_neighb_lowx)
1250 if (all(nb_ids > af_no_box))
then
1251 if (.not. any(box_done(nb_ids)) .and. &
1252 tree%boxes(nb_ids(1))%ix(1) < tree%boxes(ids(1))%ix(1) .and. &
1253 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1255 allocate(new_box_list(nx, ny))
1256 new_box_list(1, :) = nb_ids
1257 new_box_list(2:, :) = box_list
1258 box_list = new_box_list
1259 box_done(nb_ids) = .true.
1260 deallocate(new_box_list)
1265 ids = box_list(nx, :)
1266 nb_ids = tree%boxes(ids)%neighbors(af_neighb_highx)
1267 if (all(nb_ids > af_no_box))
then
1268 if (.not. any(box_done(nb_ids)) .and. &
1269 tree%boxes(nb_ids(1))%ix(1) > tree%boxes(ids(1))%ix(1) .and. &
1270 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1272 allocate(new_box_list(nx, ny))
1273 new_box_list(nx, :) = nb_ids
1274 new_box_list(1:nx-1, :) = box_list
1275 box_list = new_box_list
1276 box_done(nb_ids) = .true.
1277 deallocate(new_box_list)
1282 ids = box_list(:, 1)
1283 nb_ids = tree%boxes(ids)%neighbors(af_neighb_lowy)
1284 if (all(nb_ids > af_no_box))
then
1285 if (.not. any(box_done(nb_ids)) .and. &
1286 tree%boxes(nb_ids(1))%ix(2) < tree%boxes(ids(1))%ix(2) .and. &
1287 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1289 allocate(new_box_list(nx, ny))
1290 new_box_list(:, 1) = nb_ids
1291 new_box_list(:, 2:) = box_list
1292 box_list = new_box_list
1293 box_done(nb_ids) = .true.
1294 deallocate(new_box_list)
1299 ids = box_list(:, ny)
1300 nb_ids = tree%boxes(ids)%neighbors(af_neighb_highy)
1301 if (all(nb_ids > af_no_box))
then
1302 if (.not. any(box_done(nb_ids)) .and. &
1303 tree%boxes(nb_ids(1))%ix(2) > tree%boxes(ids(1))%ix(2) .and. &
1304 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1306 allocate(new_box_list(nx, ny))
1307 new_box_list(:, ny) = nb_ids
1308 new_box_list(:, 1:ny-1) = box_list
1309 box_list = new_box_list
1310 box_done(nb_ids) = .true.
1311 deallocate(new_box_list)
1315 if (nx == nx_prev .and. ny == ny_prev)
exit
1321 lo_bnd = af_is_phys_boundary(tree%boxes, id, af_low_neighbs)
1323 id = box_list(nx, ny)
1324 hi_bnd = af_is_phys_boundary(tree%boxes, id, af_high_neighbs)
1327 where (.not. lo_bnd) lo = lo - 1
1330 where (.not. hi_bnd) hi = hi + 1
1333 allocate(var_data(lo(1):hi(1), lo(2):hi(2), n_cc+n_add))
1337 id = box_list(ix, iy)
1339 cc(:, :, 1:n_cc) = tree%boxes(id)%cc(:, :, icc_val)
1340 if (
present(add_vars))
then
1341 call add_vars(tree%boxes(id), &
1342 cc(:, :, n_cc+1:n_cc+n_add), n_add)
1347 where ([ix, iy] == 1 .and. .not. lo_bnd) blo = 0
1350 where ([ix, iy] == [nx, ny] .and. .not. hi_bnd) bhi = nc+1
1352 vlo = blo + ([ix, iy]-1) * nc
1353 vhi = bhi + ([ix, iy]-1) * nc
1355 var_data(vlo(1):vhi(1), vlo(2):vhi(2), :) = &
1356 cc(blo(1):bhi(1), blo(2):bhi(2), :)
1361 dr = tree%boxes(id)%dr
1362 r_min = tree%boxes(id)%r_min - (1 - lo) * dr
1364 write(grid_list(i_grid),
"(A,I0)") meshdir //
'/' // grid_name, i_grid
1365 call silo_add_grid(dbix, grid_list(i_grid), 2, &
1366 hi - lo + 2, r_min, dr, 1-lo, hi - [nx, ny] * nc, n_cycle_val)
1367 write(grid_list_block(i_grid),
"(A,I0)") meshdir //
'/' // block_prefix &
1368 // grid_name, i_grid
1369 call silo_add_grid(dbix, grid_list_block(i_grid), 2, [nx+1, ny+1], &
1370 tree%boxes(id)%r_min, nc*dr, [0, 0], [0, 0], n_cycle_val)
1372 do iv = 1, n_cc+n_add
1373 write(var_list(iv, i_grid),
"(A,I0)") meshdir //
'/' // &
1374 trim(var_names(iv)) //
"_", i_grid
1375 call silo_add_var(dbix, var_list(iv, i_grid), grid_list(i_grid), &
1376 pack(var_data(:, :, iv), .true.), hi-lo+1, n_cycle_val)
1379 deallocate(var_data)
1380 deallocate(box_list)
1382 allocate(box_list(1,1,1))
1383 box_list(1,1,1) = id
1384 box_done(id) = .true.
1395 ids = pack(box_list(1, :, :), .true.)
1396 nb_ids = tree%boxes(ids)%neighbors(af_neighb_lowx)
1397 if (all(nb_ids > af_no_box))
then
1398 if (.not. any(box_done(nb_ids)) .and. &
1399 tree%boxes(nb_ids(1))%ix(1) < tree%boxes(ids(1))%ix(1) .and. &
1400 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1402 allocate(new_box_list(nx, ny, nz))
1403 new_box_list(1, :, :) = reshape(nb_ids, [ny, nz])
1404 new_box_list(2:, :, :) = box_list
1405 box_list = new_box_list
1406 box_done(nb_ids) = .true.
1407 deallocate(new_box_list)
1412 ids = pack(box_list(nx, :, :), .true.)
1413 nb_ids = tree%boxes(ids)%neighbors(af_neighb_highx)
1414 if (all(nb_ids > af_no_box))
then
1415 if (.not. any(box_done(nb_ids)) .and. &
1416 tree%boxes(nb_ids(1))%ix(1) > tree%boxes(ids(1))%ix(1) .and. &
1417 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1419 allocate(new_box_list(nx, ny, nz))
1420 new_box_list(nx, :, :) = reshape(nb_ids, [ny, nz])
1421 new_box_list(1:nx-1, :, :) = box_list
1422 box_list = new_box_list
1423 box_done(nb_ids) = .true.
1424 deallocate(new_box_list)
1429 ids = pack(box_list(:, 1, :), .true.)
1430 nb_ids = tree%boxes(ids)%neighbors(af_neighb_lowy)
1431 if (all(nb_ids > af_no_box))
then
1432 if (.not. any(box_done(nb_ids)) .and. &
1433 tree%boxes(nb_ids(1))%ix(2) < tree%boxes(ids(1))%ix(2) .and. &
1434 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1436 allocate(new_box_list(nx, ny, nz))
1437 new_box_list(:, 1, :) = reshape(nb_ids, [nx, nz])
1438 new_box_list(:, 2:, :) = box_list
1439 box_list = new_box_list
1440 box_done(nb_ids) = .true.
1441 deallocate(new_box_list)
1446 ids = pack(box_list(:, ny, :), .true.)
1447 nb_ids = tree%boxes(ids)%neighbors(af_neighb_highy)
1448 if (all(nb_ids > af_no_box))
then
1449 if (.not. any(box_done(nb_ids)) .and. &
1450 tree%boxes(nb_ids(1))%ix(2) > tree%boxes(ids(1))%ix(2) .and. &
1451 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1453 allocate(new_box_list(nx, ny, nz))
1454 new_box_list(:, ny, :) = reshape(nb_ids, [nx, nz])
1455 new_box_list(:, 1:ny-1, :) = box_list
1456 box_list = new_box_list
1457 box_done(nb_ids) = .true.
1458 deallocate(new_box_list)
1463 ids = pack(box_list(:, :, 1), .true.)
1464 nb_ids = tree%boxes(ids)%neighbors(af_neighb_lowz)
1465 if (all(nb_ids > af_no_box))
then
1466 if (.not. any(box_done(nb_ids)) .and. &
1467 tree%boxes(nb_ids(1))%ix(3) < tree%boxes(ids(1))%ix(3) .and. &
1468 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1470 allocate(new_box_list(nx, ny, nz))
1471 new_box_list(:, :, 1) = reshape(nb_ids, [nx, ny])
1472 new_box_list(:, :, 2:) = box_list
1473 box_list = new_box_list
1474 box_done(nb_ids) = .true.
1475 deallocate(new_box_list)
1480 ids = pack(box_list(:, :, nz), .true.)
1481 nb_ids = tree%boxes(ids)%neighbors(af_neighb_highz)
1482 if (all(nb_ids > af_no_box))
then
1483 if (.not. any(box_done(nb_ids)) .and. &
1484 tree%boxes(nb_ids(1))%ix(3) > tree%boxes(ids(1))%ix(3) .and. &
1485 .not. any(af_has_children(tree%boxes(nb_ids))))
then
1487 allocate(new_box_list(nx, ny, nz))
1488 new_box_list(:, :, nz) = reshape(nb_ids, [nx, ny])
1489 new_box_list(:, :, 1:nz-1) = box_list
1490 box_list = new_box_list
1491 box_done(nb_ids) = .true.
1492 deallocate(new_box_list)
1496 if (nx == nx_prev .and. ny == ny_prev .and. nz == nz_prev)
exit
1501 id = box_list(1, 1, 1)
1502 lo_bnd = af_is_phys_boundary(tree%boxes, id, af_low_neighbs)
1504 id = box_list(nx, ny, nz)
1505 hi_bnd = af_is_phys_boundary(tree%boxes, id, af_high_neighbs)
1508 where (.not. lo_bnd) lo = lo - 1
1510 hi = [nx, ny, nz] * nc
1511 where (.not. hi_bnd) hi = hi + 1
1514 allocate(var_data(lo(1):hi(1), lo(2):hi(2), lo(3):hi(3), n_cc+n_add))
1519 id = box_list(ix, iy, iz)
1520 cc(:, :, :, 1:n_cc) = tree%boxes(id)%cc(:, :, :, icc_val)
1521 if (
present(add_vars))
then
1522 call add_vars(tree%boxes(id), &
1523 cc(:, :, :, n_cc+1:n_cc+n_add), n_add)
1528 where ([ix, iy, iz] == 1 .and. .not. lo_bnd) blo = 0
1531 where ([ix, iy, iz] == [nx, ny, nz] &
1532 .and. .not. hi_bnd) bhi = nc+1
1534 vlo = blo + ([ix, iy, iz]-1) * nc
1535 vhi = bhi + ([ix, iy, iz]-1) * nc
1537 var_data(vlo(1):vhi(1), vlo(2):vhi(2), vlo(3):vhi(3), :) = &
1538 cc(blo(1):bhi(1), blo(2):bhi(2), blo(3):bhi(3), :)
1543 id = box_list(1, 1, 1)
1544 dr = tree%boxes(id)%dr
1545 r_min = tree%boxes(id)%r_min - (1 - lo) * dr
1547 write(grid_list(i_grid),
"(A,I0)") meshdir //
'/' // grid_name, i_grid
1548 call silo_add_grid(dbix, grid_list(i_grid), 3, &
1549 hi - lo + 2, r_min, dr, 1-lo, hi-[nx, ny, nz]*nc, n_cycle_val)
1550 write(grid_list_block(i_grid),
"(A,I0)") meshdir //
'/' // block_prefix // &
1552 call silo_add_grid(dbix, grid_list_block(i_grid), 3, [nx+1, ny+1, nz+1], &
1553 tree%boxes(id)%r_min, nc*dr, [0, 0, 0], [0, 0, 0], n_cycle_val)
1555 do iv = 1, n_cc+n_add
1556 write(var_list(iv, i_grid),
"(A,I0)") meshdir //
'/' // &
1557 trim(var_names(iv)) //
"_", i_grid
1558 call silo_add_var(dbix, var_list(iv, i_grid), grid_list(i_grid), &
1559 pack(var_data(:, :, :, iv), .true.), hi-lo+1, n_cycle_val)
1562 deallocate(var_data)
1563 deallocate(box_list)
1569 call silo_set_mmesh_grid(dbix, amr_name, grid_list(1:i_grid), &
1570 n_cycle_val, time_val)
1571 do iv = 1, n_cc+n_add
1572 call silo_set_mmesh_var(dbix, trim(var_names(iv)), amr_name, &
1573 var_list(iv, 1:i_grid), n_cycle_val, time_val)
1576 call silo_set_mmesh_grid(dbix, block_prefix // amr_name, &
1577 grid_list_block(1:i_grid), n_cycle_val, time_val)
1578 call silo_close_file(dbix)
1579 print *,
"af_write_silo: written " // trim(fname)
1580 end subroutine af_write_silo
1583 subroutine get_output_vars(tree, ix_out)
1584 type(af_t),
intent(in) :: tree
1585 integer,
allocatable,
intent(inout) :: ix_out(:)
1588 n = count(tree%cc_write_output(1:tree%n_var_cell))
1592 do i = 1, tree%n_var_cell
1593 if (tree%cc_write_output(i))
then
1598 end subroutine get_output_vars
This module contains the core routines of Afivo, namely those that deal with initializing and changin...
This module contains routines related to interpolation, which can interpolate 'to' the grid and 'from...
This module contains routines for writing output files with Afivo. The Silo format should probably be...
This module contains the basic types and constants that are used in the NDIM-dimensional version of A...
This file is a modification of code found in Lib_VTK_IO.
This module contains wrapper functions to simplify writing Silo files.
Type which stores all the boxes and levels, as well as some information about the number of boxes,...
The basic building block of afivo: a box with cell-centered and face centered data,...