add h5 block write

Author: fangjian

Makefile.pastr

@@ -28,7 +28,8 @@ srs=  pastr_constdef.F90 pastr_commtype.F90 pastr_commvar.F90 pastr_multiblock_t
       pastr_utility.F90 pastr_h5io.F90 pastr_io.F90 pastr_tecio.F90 pastr_xdmf.F90           \
       pastr_interpolation.F90 pastr_fsolver.F90 pastr_gradients.F90 pastr_thermo_phys.F90    \
       pastr_blasius.F90 pastr_flowvis.F90 pastr_input.F90 pastr_data_convert.F90             \
-      pastr_field_view.F90 pastr_udf.F90 pastr_process.F90 pastr_main.F90
+      pastr_field_view.F90 pastr_udf.F90 pastr_data_process.F90 pastr_process.F90            \
+      pastr_main.F90
 
 OBJS=$(srs:.F90=.o)
 

pastr/src/pastr_commtype.F90

@@ -15,11 +15,12 @@ module pastr_commtype
     end type montype
 
     type :: tblock
-      integer :: im,jm,nvar
-      integer :: ilo,ihi,jlo,jhi
+      character(len=6) :: name
+      integer :: im,jm,km,nvar
+      integer :: ilo,ihi,jlo,jhi,klo,khi
       character(len=16),allocatable :: varname(:)
-      real(wp), allocatable, dimension(:,:) :: x,y
-      real(wp), allocatable, dimension(:,:,:) :: var
+      real(wp), allocatable, dimension(:,:,:) :: x,y,z
+      real(wp), allocatable, dimension(:,:,:,:) :: var
       contains
       procedure :: init_data => alloc_block_data
     end type tblock
@@ -42,9 +43,9 @@ subroutine alloc_block_data(ablock)
       class(tblock),target :: ablock
 
       allocate(ablock%varname(ablock%nvar))
-      allocate(ablock%x(0:ablock%im,0:ablock%jm))
-      allocate(ablock%y(0:ablock%im,0:ablock%jm))
-      allocate(ablock%var(0:ablock%im,0:ablock%jm,1:ablock%nvar))
+      allocate(ablock%x(0:ablock%im,0:ablock%jm,0:ablock%km))
+      allocate(ablock%y(0:ablock%im,0:ablock%jm,0:ablock%km))
+      allocate(ablock%var(0:ablock%im,0:ablock%jm,0:ablock%km,1:ablock%nvar))
 
     end subroutine alloc_block_data
 

pastr/src/pastr_data_process.F90

@@ -0,0 +1,164 @@
+module pastr_data_process
+
+    use iso_fortran_env, only: wp => real64
+
+    implicit none
+
+    real(wp),allocatable,dimension(:,:,:) :: ro,u1,u2,u3,p,t
+    real(wp),allocatable,dimension(:,:) :: ro_xy,u1_xy,u2_xy,u3_xy,p_xy,t_xy
+
+contains
+
+    subroutine stats_read_process(mode)
+
+      use pastr_commvar, only: im,jm,km
+      use pastr_io, only : read_stats
+      
+      character(len=*),intent(in) :: mode
+
+      if(mode=='meanflow' .or. mode=='all') then
+
+        allocate( ro(0:im,0:jm,0:km),u1(0:im,0:jm,0:km),  &
+                  u2(0:im,0:jm,0:km),u3(0:im,0:jm,0:km),  &
+                  p(0:im,0:jm,0:km),t(0:im,0:jm,0:km) )
+
+        call read_stats(var=ro,varname='rom')
+        call read_stats(var=u1,varname='u1m')
+        call read_stats(var=u2,varname='u2m')
+        call read_stats(var=u3,varname='u3m')
+        call read_stats(var= p,varname= 'pm')
+        call read_stats(var= t,varname= 'tm')
+
+        allocate( ro_xy(0:im,0:jm),u1_xy(0:im,0:jm),  &
+                  u2_xy(0:im,0:jm),u3_xy(0:im,0:jm),  &
+                   p_xy(0:im,0:jm), t_xy(0:im,0:jm) )
+
+        ro_xy=average(ro)
+        u1_xy=average(u1)/ro_xy
+        u2_xy=average(u2)/ro_xy
+        u3_xy=average(u3)/ro_xy
+         p_xy=average(p)
+         t_xy=average(t)/ro_xy
+
+        call write_stats()
+
+        call plot_meanflow_xy()
+
+      elseif(mode=='2ndsta' .or. mode=='all') then
+      elseif(mode=='3rdsta' .or. mode=='all') then
+      elseif(mode=='budget' .or. mode=='all') then
+      else
+       stop 1
+      endif
+
+    end subroutine stats_read_process
+
+    subroutine write_stats
+
+      use pastr_commtype, only : tblock
+      use pastr_commvar, only: im,jm,km
+      use pastr_multiblock_type, only: block_define
+      use pastr_h5io
+
+      type(tblock),allocatable,target :: pblocks(:)
+      integer :: nblocks
+      logical :: multi_block
+      integer :: i
+      type(tblock),pointer :: b
+
+      call block_define(multi_block,nblocks,pblocks)
+
+      do i=1,nblocks
+          b=>pblocks(i)
+          b%nvar=6
+          write(b%name,'(A,I5.5)')'b',i
+          call b%init_data()
+          b%var(0:b%im,0:b%jm,0,1)=ro_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,2)=u1_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,3)=u2_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,4)=u3_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,5)= p_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,6)= t_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%varname(1) = 'ro'
+          b%varname(2) = 'u'
+          b%varname(3) = 'v'
+          b%varname(4) = 'w'
+          b%varname(5) = 'p'
+          b%varname(6) = 't'
+      enddo
+
+      call h5_write_blocks(blocks=pblocks,fname='Results/mean.fav.zm.h5')
+
+      ! call H5WriteArray(ro_xy,im,jm,'ro','Results/mean.fav.zm.h5')
+      ! call H5WriteArray(u1_xy,im,jm,'u1','Results/mean.fav.zm.h5')
+      ! call H5WriteArray(u2_xy,im,jm,'u2','Results/mean.fav.zm.h5')
+      ! call H5WriteArray(u3_xy,im,jm,'u3','Results/mean.fav.zm.h5')
+      ! call H5WriteArray( p_xy,im,jm, 'p','Results/mean.fav.zm.h5')
+      ! call H5WriteArray( t_xy,im,jm, 't','Results/mean.fav.zm.h5')
+
+    end subroutine write_stats
+
+    subroutine plot_meanflow_xy
+
+      use pastr_commvar, only: im,jm,km
+      use pastr_io, only : read_grid
+      use pastr_multiblock_type, only: block_define
+      use pastr_tecio
+      use pastr_commtype, only : tblock
+
+      real(wp),allocatable :: x(:,:),y(:,:),z(:,:)
+      type(tblock),allocatable,target :: pblocks(:)
+      integer :: nblocks
+      logical :: multi_block
+      integer :: i
+      type(tblock),pointer :: b
+
+      allocate(x(0:im,0:jm),y(0:im,0:jm),z(0:im,0:jm))
+
+      call read_grid(x=x,y=y,z=z,kslice=0)
+
+      call block_define(multi_block,nblocks,pblocks)
+
+      do i=1,nblocks
+          b=>pblocks(i)
+          b%nvar=6
+          call b%init_data()
+          b%x(0:b%im,0:b%jm,0)=x(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%y(0:b%im,0:b%jm,0)=y(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,1)=ro_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,2)=u1_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,3)=u2_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,4)=u3_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,5)= p_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,6)= t_xy(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%varname(1) = 'ro'
+          b%varname(2) = 'u'
+          b%varname(3) = 'v'
+          b%varname(4) = 'w'
+          b%varname(5) = 'p'
+          b%varname(6) = 't'
+
+      enddo
+
+      call tecbin(filename='Results/tecmean',block=pblocks)
+
+    end subroutine plot_meanflow_xy
+
+    function average(var) result(varm)
+
+      use pastr_commvar, only: im,jm,km
+
+      real(wp),intent(in) :: var(0:im,0:jm,0:km)
+      real(wp) :: varm(0:im,0:jm)
+
+      integer :: k
+
+      varm=0._wp
+      do k=1,km
+        varm(:,:)=varm(:,:)+var(:,:,k)
+      enddo
+      varm=varm/dble(km)
+
+    end function average
+
+end module pastr_data_process

pastr/src/pastr_field_view.F90

@@ -14,6 +14,7 @@ module pastr_field_view
     subroutine write_xy_slice(filein,fileout,nfirst,nlast,slice,format)
 
       use pastr_io, only : read_grid,parse_command_line
+      use pastr_multiblock_type, only: block_define
       use pastr_commvar, only : im,jm,km
       use pastr_h5io
       use pastr_tecio
@@ -101,9 +102,9 @@ subroutine write_xy_slice(filein,fileout,nfirst,nlast,slice,format)
           b=>pblocks(i)
           b%nvar=num_var_out
           call b%init_data()
-          b%x(0:b%im,0:b%jm)=x(b%ilo:b%ihi,b%jlo:b%jhi)
-          b%y(0:b%im,0:b%jm)=y(b%ilo:b%ihi,b%jlo:b%jhi)
-          b%var(0:b%im,0:b%jm,1:num_var_out)=dat_var_out(b%ilo:b%ihi,b%jlo:b%jhi,1:num_var_out )
+          b%x(0:b%im,0:b%jm,0)=x(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%y(0:b%im,0:b%jm,0)=y(b%ilo:b%ihi,b%jlo:b%jhi)
+          b%var(0:b%im,0:b%jm,0,1:num_var_out)=dat_var_out(b%ilo:b%ihi,b%jlo:b%jhi,1:num_var_out )
           b%varname=nam_var_out
         enddo
 
@@ -387,40 +388,4 @@ subroutine var_define(num_var_in,nam_var_in,num_var_out,nam_var_out)
 
     end subroutine var_define
 
-    subroutine block_define(multi_block,nblocks,pblocks)
-
-      use pastr_commtype
-
-      logical,intent(out) :: multi_block
-      integer,intent(out) :: nblocks
-      type(tblock),intent(out),allocatable :: pblocks(:)
-
-      integer :: i
-      logical :: lfex
-
-      inquire(file='blockdef.txt',exist=lfex)
-      if(lfex) then
-        multi_block=.true.
-      else
-        multi_block=.false.
-        return
-      endif
-
-      open(12,file='blockdef.txt')
-      read(12,*)nblocks
-      allocate(pblocks(nblocks))
-      do i=1,nblocks
-        read(12,*)pblocks(i)%ilo,pblocks(i)%ihi,pblocks(i)%jlo,pblocks(i)%jhi
-        pblocks(i)%im=pblocks(i)%ihi-pblocks(i)%ilo
-        pblocks(i)%jm=pblocks(i)%jhi-pblocks(i)%jlo
-      enddo
-      close(12)
-      print*,' >> blockdef.txt'
-
-      print*,' ** nblocks: ',nblocks
-
-      return
-
-    end subroutine block_define
-
 end module pastr_field_view

pastr/src/pastr_fsolver.F90

@@ -1427,6 +1427,7 @@ subroutine fsolve ( fcn, n, x, fvec, tol, info )
     !*****************************************************************************80
     !
     !! fsolve() seeks a zero of N nonlinear equations in N variables.
+    !! Return the roots of the (non-linear) equations defined by fcn(x) = 0 given a starting estimate.
     !
     !  Discussion:
     !

pastr/src/pastr_h5io.F90

@@ -1650,7 +1650,59 @@ end subroutine h5_writesubset
   !+-------------------------------------------------------------------+
   ! this subroutine is used to read a 3d array file with hdf5
   !+-------------------------------------------------------------------+
-  
+
+  subroutine h5_write_blocks(blocks,fname)
+    
+    use pastr_commtype, only : tblock
+
+    type(tblock),intent(in),target :: blocks(:)
+    character(len=*),intent(in) :: fname
+
+    logical :: lfilalive
+    integer :: i,n
+    type(tblock),pointer :: b
+    integer(hid_t) :: file_id,group_id
+    ! file identifier
+    integer(hid_t) :: dataspace_id,dset_id
+    ! dataset identifier
+    integer :: h5error ! error flag
+    integer(hsize_t) :: dimt(3)
+    
+    call h5open_f(h5error)
+    !
+    inquire(file=fname, exist=lfilalive)
+    if(lfilalive) then
+      call h5fopen_f(fname,H5F_ACC_RDWR_F,file_id,h5error)
+    else
+      call h5fcreate_f(fname,H5F_ACC_TRUNC_F,file_id,h5error)
+    end if
+    if(h5error.ne.0)  stop ' !! error in h5_write_blocks 1'
+    
+    do i=1,size(blocks)
+
+      b=>blocks(i)
+
+      call h5gcreate_f(file_id, b%name, group_id, h5error)
+
+      dimt=(/b%im+1,b%jm+1,b%km+1/)
+
+      do n=1,b%nvar
+        CALL h5screate_simple_f(3, dimt, dataspace_id, h5error)
+        call h5dcreate_f(group_id, trim(b%varname(n)), h5t_native_double, &
+                         dataspace_id, dset_id, h5error)
+        call h5dwrite_f(dset_id, h5t_native_double, b%var(:,:,:,n), dimt, h5error)
+        CALL h5sclose_f(dataspace_id, h5error)
+        call h5dclose_f(dset_id, h5error)
+      enddo
+      call h5gclose_f(group_id, h5error)
+
+    enddo
+    call h5fclose_f(file_id, h5error)
+
+    call h5close_f(h5error)
+
+  end subroutine h5_write_blocks
+
   subroutine h5_writearray3d(varin,dim1,dim2,dim3,vname,fname,explicit,ierr)
     !
     integer,intent(in) :: dim1,dim2,dim3

pastr/src/pastr_io.F90

@@ -8,6 +8,7 @@ module pastr_io
 
     public :: parse_command_line,read_grid
     public :: global_n_arguments,read_monitor_data,write_monitor_data
+    public :: read_stats
 
     integer :: global_n_arguments=0
 
@@ -196,4 +197,27 @@ subroutine write_monitor_data(mon_data)
 
     end subroutine write_monitor_data
 
+    subroutine read_stats(var,varname)
+
+      use pastr_commvar, only: im,jm,km
+      use pastr_h5io
+
+      real(wp) :: var(0:im,0:jm,0:km)
+      character(len=*),intent(in) :: varname
+
+      character(len=18) :: fname
+      integer,save :: nsamples=0
+
+      fname='outdat/meanflow.h5'
+
+      if(nsamples==0) then
+        call H5ReadArray(nsamples,'nsamples',fname)
+        print*,' ** nsamples: ',nsamples
+      endif
+
+      call H5ReadArray(var,im,jm,km,varname,fname)
+      var = var/dble(nsamples)
+
+    end subroutine read_stats
+
 end module pastr_io