add xdmf output to help loading hdf5 file in paraview

Makefile

@@ -95,7 +95,10 @@ $(OBJDIR)/%.o: %.f90 $(MOBJS)
 	$(FC) $(FFLAGS) -c -o $@ $< $(LDFLAGS)
 
 clean:
-	rm -rf $(PROGRAMS) $(OBJDIR)/*.o $(OBJDIR)/*.mod a.out
+	rm -rf $(PROGRAMS) $(OBJDIR)/*.o $(OBJDIR)/*.mod UP2D
+
+cleanall:
+	\rm -rf *.t *.h5 *.xmf
 
 tidy:
 	rm -rf $(OBJDIR)/*.o $(OBJDIR)/*.mod a.out

src/fileIO/save_fields.f90

@@ -1,5 +1,145 @@
+module utils
+
+  implicit none
+
+contains
+
+  function make_timestring(time) result(res)
+
+  use vars
+  implicit none
+
+  character(:), allocatable  :: res
+  real(kind=pr), intent (in) :: time
+
+  character(len=strlen)      :: tmp
+  write(tmp,'(i6.6)') nint(time*100.d0)
+
+  res = trim(tmp)
+
+end function make_timestring
+
+end module utils
+
+subroutine write_xdmf(time, it)
+  use vars
+  use utils
+  implicit none
+
+  real(kind=pr), intent (in) :: time
+  integer,       intent (in) :: it
+
+  character(LEN=80) :: filename
+  character(LEN=80) :: filebase
+  character(LEN=80) :: charbuf, char_nx_ny
+  character(1), parameter :: endl  = char(10)  ! end of line
+
+  integer           :: error
+
+  integer, parameter    :: IO_NODES_PER_CELL=4
+  character(len=80), parameter :: IO_TOPOLOGY_TYPE = 'Quadrilateral'
+
+  filebase = 'UP2D_' // make_timestring(time)
+  filename = 'UP2D_' // make_timestring(time) // '.xmf'
+
+  open(10,file=filename,iostat=error)
+
+  write(10,'(a)') '<?xml version="1.0" ?>'
+  write(10,'(a)') '<Xdmf xmlns:xi="http://www.w3.org/2003/XInclude" Version="2.2">'
+  write(10,'(a)') '  <Domain>'
+  write(10,'(a)') '    <Grid Name="' // trim(filebase) // '" GridType="Uniform">'
+
+  write(charbuf,fmt='(1(f7.2))',iostat=error) time
+  write(10,'(a)') '      <Time TimeType="Single" Value="' // trim(charbuf) // '" />'
+
+  ! Topology
+  write(char_nx_ny,fmt='((i7,a,i7))',iostat=error) nx,' ',ny
+  write(10,'(a)') '      <Topology TopologyType="2DCoRectMesh" NumberOfElements="' // trim(char_nx_ny) // '"  />'
+
+  write(10,'(a)') '    <Geometry Type="ORIGIN_DXDY">'
+  write(10,'(a)') '    <DataStructure'
+  write(10,'(a)') '       Name="Origin"'
+  write(10,'(a)') '       DataType="' // 'Double'  // '"'
+  write(10,'(a)') '       Dimensions="2"'
+  write(10,'(a)') '       Format="XML">'
+  write(10,'(a)') '       0 0'
+  write(10,'(a)') '    </DataStructure>'
+  write(10,'(a)') '    <DataStructure'
+  write(10,'(a)') '       Name="Spacing"'
+  write(10,'(a)') '       DataType="' // 'Double' // '"'
+  write(10,'(a)') '       Dimensions="2"'
+  write(10,'(a)') '       Format="XML">'
+  write(10,'(a)') '       1 1'
+  write(10,'(a)') '    </DataStructure>'
+  write(10,'(a)') '    </Geometry>' 
+
+  ! write data attributes
+  if ( iSaveVorticity == 1) then
+     write(10,'(a)') '      <Attribute Center="Node" Name="' // 'vor' // '">'
+     write(10,'(a)') '        <DataStructure'
+     write(10,'(a)') '           DataType="' // 'Double' // '"'
+     write(10,'(a)') '           Dimensions="' // trim(char_nx_ny) // '"'
+     write(10,'(a)') '           Format="HDF">'
+     write(10,'(a)') '           ' // 'vor_' // make_timestring(time) // '.h5' // ':/vor'
+     write(10,'(a)') '        </DataStructure>'
+     write(10,'(a)') '      </Attribute>'
+  end if
+
+  if ( iSaveVelocity == 1) then
+     write(10,'(a)') '      <Attribute Center="Node" Name="' // 'ux' // '">'
+     write(10,'(a)') '        <DataStructure'
+     write(10,'(a)') '           DataType="' // 'Double' // '"'
+     write(10,'(a)') '           Dimensions="' // trim(char_nx_ny) // '"'
+     write(10,'(a)') '           Format="HDF">'
+     write(10,'(a)') '           ' // 'ux_' // make_timestring(time) // '.h5' // ':/ux'
+     write(10,'(a)') '        </DataStructure>'
+     write(10,'(a)') '      </Attribute>'
+
+     write(10,'(a)') '      <Attribute Center="Node" Name="' // 'uy' // '">'
+     write(10,'(a)') '        <DataStructure'
+     write(10,'(a)') '           DataType="' // 'Double' // '"'
+     write(10,'(a)') '           Dimensions="' // trim(char_nx_ny) // '"'
+     write(10,'(a)') '           Format="HDF">'
+     write(10,'(a)') '           ' // 'uy_' // make_timestring(time) // '.h5' // ':/uy'
+     write(10,'(a)') '        </DataStructure>'
+     write(10,'(a)') '      </Attribute>'
+  end if
+
+  if ( iSaveMask == 1) then
+     write(10,'(a)') '      <Attribute Center="Node" Name="' // 'mask' // '">'
+     write(10,'(a)') '        <DataStructure'
+     write(10,'(a)') '           DataType="' // 'Double' // '"'
+     write(10,'(a)') '           Dimensions="' // trim(char_nx_ny) // '"'
+     write(10,'(a)') '           Format="HDF">'
+     write(10,'(a)') '           ' // 'mask_' // make_timestring(time) // '.h5' // ':/mask'
+     write(10,'(a)') '        </DataStructure>'
+     write(10,'(a)') '      </Attribute>'
+  end if
+
+  if ( iSavePressure == 1 ) then
+     write(10,'(a)') '      <Attribute Center="Node" Name="' // 'p' // '">'
+     write(10,'(a)') '        <DataStructure'
+     write(10,'(a)') '           DataType="' // 'Double' // '"'
+     write(10,'(a)') '           Dimensions="' // trim(char_nx_ny) // '"'
+     write(10,'(a)') '           Format="HDF">'
+     write(10,'(a)') '           ' // 'p_' // make_timestring(time) // '.h5' // ':/p'
+     write(10,'(a)') '        </DataStructure>'
+     write(10,'(a)') '      </Attribute>'
+  end if
+
+  ! write footer
+  write(10,'(a)') '    </Grid>'
+  write(10,'(a)') '  </Domain>'
+  write(10,'(a)') '</Xdmf>'
+
+
+  close(10)
+
+end subroutine write_xdmf
+
 subroutine save_fields(time, it, u, uk, vort, mask, us, mask_sponge)
   use vars
+  use utils
   implicit none
 
   real(kind=pr), intent (in) :: time
@@ -10,12 +150,14 @@ subroutine save_fields(time, it, u, uk, vort, mask, us, mask_sponge)
   real(kind=pr),dimension(0:nx-1,0:ny-1,1:2),intent(inout) :: us
 
   real(kind=pr),dimension(:,:), allocatable :: work, pk
-  character(len=strlen) :: timestring
+  character(:),allocatable :: timestring
 
   allocate(work(0:nx-1, 0:ny-1), pk(0:nx-1, 0:ny-1))
 
-  write(timestring,'(i6.6)') nint(time*100.d0)
+  !write(timestring,'(i6.6)') nint(time*100.d0)
   !write(timestring,'(i6.6)') it
+  timestring=make_timestring(time)
+
   write(*,'("Saving. time=",es12.4," vormax=",es12.4," fname=",A)') time, maxval(vort), timestring
 
   if ( iSaveVorticity == 1) then

src/time_step.f90

@@ -41,7 +41,8 @@ subroutine time_step (u, uk, nlk, pk, vort, mask, us, mask_sponge)
 
     if ( time_for_output( time, dt1, it, tsave, itsave, Tmax, 0.d0) ) then
       ! save output fields to disk
-      call save_fields(time, it, u, uk, vort, mask, us, mask_sponge)
+       call save_fields(time, it, u, uk, vort, mask, us, mask_sponge)
+       call write_xdmf(time, it)
     endif
 
     ! output remaining time