subroutine write_vtu_binary(params, flow, mesh, fields, species, energy, transport, step, time)
type(case_params_t), intent(in) :: params
type(flow_mpi_t), intent(inout) :: flow
type(mesh_t), intent(in) :: mesh
type(flow_fields_t), intent(in) :: fields
type(species_fields_t), intent(in) :: species
type(energy_fields_t), intent(in) :: energy
type(transport_properties_t), intent(in) :: transport
integer, intent(in) :: step
real(rk), intent(in) :: time
integer :: unit_id
integer :: p, c, n, m, k, n_owned
character(len=path_len + 64) :: filename, local_filename
integer, allocatable :: owned_indices(:)
integer(8) :: cur_offset
integer(8) :: off_vel, off_umag, off_pres, off_thermo, off_div, off_rho, off_nu
integer(8), allocatable :: off_species(:)
integer(8) :: off_cp, off_lambda, off_temp, off_enth, off_qrad
integer(8) :: off_points, off_conn, off_offsets, off_types
character(len=2048) :: line_buf
real(rk), allocatable :: r_buf(:)
real(rk), allocatable :: velocity_buf(:,:)
real(rk), allocatable :: points_buf(:,:)
integer(4), allocatable :: conn_buf(:,:)
integer(4), allocatable :: off_buf(:)
integer(4), allocatable :: type_buf(:)
! Count owned cells
n_owned = 0
do c = 1, mesh%ncells
if (flow%cell_owner(c) == flow%rank) n_owned = n_owned + 1
end do
allocate(owned_indices(n_owned))
n_owned = 0
do c = 1, mesh%ncells
if (flow%cell_owner(c) == flow%rank) then
n_owned = n_owned + 1
owned_indices(n_owned) = c
end if
end do
! Compute Offsets
cur_offset = 0
off_vel = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*3*8
off_umag = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
off_pres = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
off_thermo = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
off_div = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
off_rho = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
off_nu = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
if (params%enable_species .and. species%nspecies > 0) then
allocate(off_species(species%nspecies))
do k = 1, species%nspecies
off_species(k) = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
end do
end if
if (params%enable_energy) then
off_cp = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
off_lambda = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
off_temp = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
off_enth = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
off_qrad = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8
end if
off_points = cur_offset
cur_offset = cur_offset + 8 + int(mesh%npoints, 8)*3*8
off_conn = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*8*4
off_offsets = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*4
off_types = cur_offset
cur_offset = cur_offset + 8 + int(n_owned, 8)*4
! Open stream binary file
write(local_filename,'("flow_",i6.6,"_P",i4.4,".vtu")') step, flow%rank
filename = trim(params%output_dir)//'/VTK/'//trim(local_filename)
open(newunit=unit_id, file=trim(filename), access='stream', form='unformatted', status='replace')
! Write XML header and declarations
write(unit_id) '<?xml version="1.0"?>' // char(10)
write(unit_id) '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">' // char(10)
write(unit_id) ' <UnstructuredGrid>' // char(10)
write(line_buf, '(" <Piece NumberOfPoints=""",i0,""" NumberOfCells=""",i0,""">")') mesh%npoints, n_owned
write(unit_id) trim(line_buf) // char(10)
write(unit_id) ' <PointData>' // char(10)
write(unit_id) ' </PointData>' // char(10)
write(unit_id) ' <CellData Scalars="pressure" Vectors="velocity">' // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""velocity"" NumberOfComponents=""3"" format=""appended"" offset=""",i0,"""/>")') off_vel
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""U_mag"" format=""appended"" offset=""",i0,"""/>")') off_umag
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""pressure"" format=""appended"" offset=""",i0,"""/>")') off_pres
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""thermo_pressure"" format=""appended"" offset=""",i0,"""/>")') off_thermo
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""divergence"" format=""appended"" offset=""",i0,"""/>")') off_div
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""rho"" format=""appended"" offset=""",i0,"""/>")') off_rho
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""nu"" format=""appended"" offset=""",i0,"""/>")') off_nu
write(unit_id) trim(line_buf) // char(10)
if (params%enable_species .and. species%nspecies > 0) then
do k = 1, species%nspecies
write(line_buf, '(" <DataArray type=""Float64"" Name=""Y_",a,""" format=""appended"" offset=""",i0,"""/>")') &
trim(params%species_name(k)), off_species(k)
write(unit_id) trim(line_buf) // char(10)
end do
end if
if (params%enable_energy) then
write(line_buf, '(" <DataArray type=""Float64"" Name=""cp"" format=""appended"" offset=""",i0,"""/>")') off_cp
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""lambda"" format=""appended"" offset=""",i0,"""/>")') off_lambda
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""temperature"" format=""appended"" offset=""",i0,"""/>")') off_temp
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""enthalpy"" format=""appended"" offset=""",i0,"""/>")') off_enth
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""qrad"" format=""appended"" offset=""",i0,"""/>")') off_qrad
write(unit_id) trim(line_buf) // char(10)
end if
write(unit_id) ' </CellData>' // char(10)
write(unit_id) ' <Points>' // char(10)
write(line_buf, '(" <DataArray type=""Float64"" Name=""Points"" NumberOfComponents=""3"" format=""appended"" offset=""",i0,"""/>")') off_points
write(unit_id) trim(line_buf) // char(10)
write(unit_id) ' </Points>' // char(10)
write(unit_id) ' <Cells>' // char(10)
write(line_buf, '(" <DataArray type=""Int32"" Name=""connectivity"" format=""appended"" offset=""",i0,"""/>")') off_conn
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Int32"" Name=""offsets"" format=""appended"" offset=""",i0,"""/>")') off_offsets
write(unit_id) trim(line_buf) // char(10)
write(line_buf, '(" <DataArray type=""Int32"" Name=""types"" format=""appended"" offset=""",i0,"""/>")') off_types
write(unit_id) trim(line_buf) // char(10)
write(unit_id) ' </Cells>' // char(10)
write(unit_id) ' </Piece>' // char(10)
write(unit_id) ' </UnstructuredGrid>' // char(10)
! Write the AppendedData raw block header
write(unit_id) ' <AppendedData encoding="raw">' // char(10)
write(unit_id) '_'
! Allocate reusable cell buffer
allocate(r_buf(n_owned))
! 1. velocity
allocate(velocity_buf(3, n_owned))
do n = 1, n_owned
c = owned_indices(n)
velocity_buf(:,n) = fields%u(:,c)
end do
write(unit_id) int(n_owned * 3 * 8, 8)
write(unit_id) velocity_buf
deallocate(velocity_buf)
! 2. U_mag
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = sqrt(dot_product(fields%u(:,c), fields%u(:,c)))
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
! 3. pressure
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = fields%p(c)
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
! 4. thermo_pressure
r_buf = params%background_press
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
! 5. divergence
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = fields%div(c)
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
! 6. rho
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = transport%rho(c)
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
! 7. nu
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = transport%nu(c)
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
! 8. species
if (params%enable_species .and. species%nspecies > 0) then
do k = 1, species%nspecies
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = species%Y(k,c)
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
end do
deallocate(off_species)
end if
! 9. cp, lambda, temp, enth, qrad (if energy)
if (params%enable_energy) then
! cp
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = energy%cp(c)
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
! lambda
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = energy%lambda(c)
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
! temperature
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = energy%T(c)
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
! enthalpy
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = energy%h(c)
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
! qrad
do n = 1, n_owned
c = owned_indices(n)
r_buf(n) = energy%qrad(c)
end do
write(unit_id) int(n_owned * 8, 8)
write(unit_id) r_buf
end if
deallocate(r_buf)
! 10. Points
allocate(points_buf(3, mesh%npoints))
do n = 1, mesh%npoints
points_buf(:,n) = mesh%points(:,n)
end do
write(unit_id) int(mesh%npoints * 3 * 8, 8)
write(unit_id) points_buf
deallocate(points_buf)
! 11. Connectivity
allocate(conn_buf(8, n_owned))
do n = 1, n_owned
c = owned_indices(n)
do k = 1, 8
conn_buf(k, n) = int(mesh%cells(c)%nodes(k) - 1, 4)
end do
end do
write(unit_id) int(n_owned * 8 * 4, 8)
write(unit_id) conn_buf
deallocate(conn_buf)
! 12. Offsets
allocate(off_buf(n_owned))
do n = 1, n_owned
off_buf(n) = int(n * 8, 4)
end do
write(unit_id) int(n_owned * 4, 8)
write(unit_id) off_buf
deallocate(off_buf)
! 13. Types
allocate(type_buf(n_owned))
type_buf = 12_4
write(unit_id) int(n_owned * 4, 8)
write(unit_id) type_buf
deallocate(type_buf)
! Close the raw block tag and file
write(unit_id) char(10) // ' </AppendedData>' // char(10)
write(unit_id) '</VTKFile>' // char(10)
close(unit_id)
deallocate(owned_indices)
end subroutine write_vtu_binary