Hello,
I am trying to visualize 3D arrays over time that I am generating in a Fortran code. Do you know of any publicly accessible libraries or code that can write 3D arrays in a format that can be view using VisIt.
1 Like
Hi @Woolford180 welcome to the Discourse.
Haven’t used VisIt but read about a few months ago and saw it supports VTK, so you could try out GitHub - szaghi/VTKFortran: pure Fortran VTK (XML) API
@Woolford180, welcome to the forum!
In the past I have also written the text based VTK format from Fortran (hfsolver/src/fe_mesh.f90 at b4c50c1979fb7e468b1852b144ba756f5a51788d · certik/hfsolver · GitHub), and you can then convert it to more modern VTK format with a script, and then load into Paraview or Visit.
I use VisIt as my main postprocessor for my CFD code, and I use nowadays VTK exclusively:
https://bazaar.launchpad.net/~hwkrus/dolfyn-cfd/trunk/view/head:/src/vtk.f90
This was done using VTK and VisIt (the ground being a Google Maps picture):
2 Likes
Just be aware of which of the VTK formats you are going to use (text is always the safest if you don’t have super large files). If you use the old “legacy” format and want to write out data as binary instead of text you HAVE to write the data in BIG endian format. Fortunately most compilers will let you specify conversion to BIG endian as an option on the OPEN statement. I don’t know about the newer XML based format. You can embed binary into the XML but it has to be base64 encoded if written by Fortran. I think C can embed raw binary (as a stream) but I’m not sure how you do that in Fortran.
Edit. Here are the open statements I use in my Legacy VTK writer package. Note this works with ifort since thats really all I use but should work with other compilers that support CONVERT type keywords on OPEN statements
If (fileType == VTK_ASCII) Then
Open(NEWUNIT=vtkunit, FILE=thisFile, ACCESS="SEQUENTIAL", &
FORM="FORMATTED", STATUS="UNKNOWN")
Else !must have convert to big endian for binary to work
If (check_endian() == LITTLE_ENDIAN) Then
If (print_VTK_warnings) Then
Print *,''
Print *,' VTKwriteOpenFile - WARNING:'
Print *,''
Print *,' This processor uses LITTLE_ENDIAN byte ordering.'
Print *,' Legacy VTK requires BIG_ENDIAN ordering. Therefore,'
Print *,' this file will be created with the CONVERT attribute'
Print *,' set to BIG_ENDIAN'
Print *,''
End If
Open(NEWUNIT=vtkunit, FILE=thisFile, ACCESS="STREAM", &
FORM="UNFORMATTED", ACTION="WRITE", CONVERT="BIG_ENDIAN", &
STATUS="UNKNOWN")
Else
Open(NEWUNIT=vtkunit, FILE=thisFile, ACCESS="STREAM", &
FORM="UNFORMATTED", ACTION="WRITE", STATUS="UNKNOWN")
EndIf
accessType = "STREAM"
EndIf
Note you will need a utility to check the endianess of you current computer.
1 Like
@Hendrik it looks like your code also uses the text based VTK. I had good experience with that also.
try this
program vtk_3d
implicit none
integer , parameter :: Nx = 64
integer , parameter :: Ny = 64
integer , parameter :: Nz = 64
integer , parameter :: dx = 1
integer , parameter :: dy = 1
integer , parameter :: dz = 1
integer , parameter :: nsteps = 1000
integer :: tsteps
real , dimension ( Nx, Ny, Nz ) :: r, var
integer :: i, j, k
!iterations
integer, parameter :: frequency = 100
character ( len = 80 ) :: filename
integer :: lun
integer :: npoint
real :: x, y, z
! start loop
time_loop:do tsteps = 1, nsteps
call random_number ( r )
do i = 1, Nx
do j = 1, Ny
do k = 1, Nz
var(i,j,k) = var(i,j,k) + r(i,j,j)
end do
end do
end do
! saving vtk files
npoint = Nx*Ny*Nz
vtk:if ( mod ( tsteps, frequency ).eq.0 ) then
write ( filename, '( "test_", i0.3, ".vtk" )') tsteps
open ( newunit = lun, file = filename, status = 'replace' )
! header of VTK file
write ( lun, '(a)' ) '# vtk DataFile Version 2.0'
write ( lun, '(a)' ) 'time_test.vtk'
write ( lun, '(a)' ) 'ASCII'
write ( lun, '(a)' ) 'DATASET STRUCTURED_GRID'
! coordinates of grid points
write ( lun, '(a, 3i6)' ) 'DIMENSIONS', Nx, Ny, Nz
write ( lun, '(ai6, 2x, a)' ) 'POINTS ', npoint, 'float'
do i = 1, Nx
do j = 1, Ny
do k = 1, Nz
x = ( i - 1 )*dx
y = ( j - 1 )*dy
z = ( k - 1 )*dz
write ( lun, '(2x, g14.6, 2x, g14.6, 2x, g14.6)' ) x, y, z
end do
end do
end do
! write grid point values
write ( lun, '(a,i6)' ) 'POINT_DATA ', NPOINT
write ( lun, '(a)' ) 'SCALARS var float'
write ( lun, '(a)' ) 'LOOKUP_TABLE default'
do i = 1, Nx
do j = 1, Ny
do k = 1, Nz
write ( lun, '(1x, g14.6)' ) var(i,j,k)
end do
end do
end do
close ( unit = lun )
endif vtk
end do time_loop
end program vtk_3d
Yes it does. It is simpler. I am not a fan of xml.
It caters for big and little endian, ascii or binary (ascii is mainly for testing purposes). So it relies on Fortran 2003 and up (streaming files). Used it for models containing many millions of cells, nodes and data. (You’ll find some commented out testing hacks in the code though).
In the picture you’ll see a yellow line, a glidepath or approach. Things like that (or a runway) I simply write by hand:
# vtk DataFile Version 2.0
Approach
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 2 float
-978.96 72.63 86.38 403.43 970.37 0.0
CELLS 1 3
2 0 1
CELL_TYPES 1
3
POINT_DATA 2
SCALARS kleur int
LOOKUP_TABLE default
1 1
I really like VisIt, as you can write scripts for it in Python.
1 Like
These guys (MPI-AMRVAC - Astrophysical MHD simulations) use vtk as one of the potential output format. I normally these and analyze with VisIt afterwards. I use mainly binary vtk files. Which can be, in addition, MPI-oriented, i.e. bigger files are split in pieces such that reading and processing by VisIt can be done in parallel. Note for that you will need to recompile VisIt with parallel support also.