Skip to content

Commit

Permalink
Routine for computing 1dspectra
Browse files Browse the repository at this point in the history
  • Loading branch information
Kevin Zhong committed Dec 13, 2024
1 parent ca9afc0 commit d50199e
Show file tree
Hide file tree
Showing 10 changed files with 219 additions and 20 deletions.
2 changes: 1 addition & 1 deletion HITForcing.in
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
WHICH_HIT
1
TL EPSSTAR KF / KMIN
0.399 0.05 2.3
0.014 0.01 2.3
C
0.1
8 changes: 4 additions & 4 deletions bou.in
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,12 @@ XLEN YLEN ZLEN
NREAD PREAD
0 0
NTST NSST TFRAME TPIN IRESET
6000 3 1e-2 5e-3 0
6000 3 1e-1 5e-3 0
IDTV DT DTMAX WANT_CFL
0 1e-2 1e-2 1.0
1 1e-2 1e-2 1.0
REN PRANDTL BETA_GZ
50.0 1.0 0.0
650.0 1.0 0.0
TMELT TLIQUID TSOLID LATHEAT CP_LIQUID TEMP_RESTART
0.0d0 1.0d0 0.0d0 4.0d0 1.0d0 0
IHIT
0
1
14 changes: 9 additions & 5 deletions gcurv.f90
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@ subroutine gcurv
call get_prow_pcol ! KZ: pencils for bounding box in ray-tagging
call InitArrays

if (specflag) call initSpectra

tin(1) = MPI_WTIME()

call HdfStart
Expand Down Expand Up @@ -55,8 +57,8 @@ subroutine gcurv
cflm=0.d0

!call inqpr_rotated
call inqpr
!call inqpr_taylorGreen
!call inqpr
call inqpr_taylorGreen

else

Expand Down Expand Up @@ -157,8 +159,10 @@ subroutine gcurv
call mkmov_hdf_zcut
call write_tecplot_geom
!call mpi_write_tempField
!call mpi_write_vel
call mpi_write_vel
!call mpi_write_field

if (specflag) call compute_1d_spectra
endif

! ASCII write
Expand All @@ -176,8 +180,8 @@ subroutine gcurv
call CalcDissipation

! KZ: relative Lagrangian motion tracking
call calcRelBoxVel
call calcRelShellVel
!call calcRelBoxVel
!call calcRelShellVel

time=time+dt
!if((ntime.eq.ntst).or.(mod(ntime,1000).eq.0)) then !to perform when needed not only at the end
Expand Down
7 changes: 7 additions & 0 deletions init.f90
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,13 @@ subroutine InitArrays
call AllocateReal3DArray(d_UsolidT_dxj,1,n1,1,n2,kstart,kend)
call AllocateLogical3DArray(solid_mask,1,n1,1,n2,kstart,kend)

! Default values for single-phase
VOFx(:,:,:) = 1.
VOFy(:,:,:) = 1.
VOFz(:,:,:) = 1.
VOFp(:,:,:) = 1.
solid_mask(:,:,:) = .false.

! Auxilary fractional-step pseudo-pressure
call AllocateReal3DArray(dph,1,n1,1,n2+1, &
kstart-lvlhalo,kend+lvlhalo)
Expand Down
13 changes: 8 additions & 5 deletions makefile
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,18 @@ ifdef TIMED
FC += -DTIMED
endif

# Flag for computing spectra, do "make SPEC=1"
ifdef SPEC
FC += -DSPEC
endif

FFILES = auxroutines.f90 cfl.f90 cordin.f90 divg.f90 gcurv.f90 hdf.f90 \
hdnl1.f90 hdnl2.f90 hdnl3.f90 hdnlte.f90 hit.f90 inirea.f90 init.f90 inqpr.f90 \
interp.f90 invtr1.f90 invtr2.f90 invtr3.f90 invtrte.f90 matrix_transpose.f90 \
mpi_routines.f90 mpiauxroutines.f90 papero.f90 phcalc.f90 phini.f90 \
prcalc.f90 quit.f90 solxi.f90 solxj.f90 solxk.f90 stat.f90 \
tridiag_periodic.f90 tsch.f90 updvp.f90 inicut.f90 movcut.f90 hdf2.f90 \
diss.f90 vorticity.f90 injection.f90 calcSlipVels.f90 \
diss.f90 vorticity.f90 injection.f90 calcSlipVels.f90 spectra.f90 \
sphereTagging.f90

FFILES += allotri.f90 RigidAuxRoutines.f90 create_geo.f90 findCentroidIndices.f90 remesh_coarsen.f90 remesh_smooth.f90 findProbeIndices.f90 \
Expand Down Expand Up @@ -66,9 +71,7 @@ $(OBJDIR)/%.o: %.f90 $(MOBJS)


#-- Output folders
OUTDIR := flowmov continuation stringdata


OUTDIR := flowmov continuation stringdata spectra

clean :
rm -rf $(OBJDIR)
Expand All @@ -88,4 +91,4 @@ ${OBJDIR}:

outdir: ${OUTDIR}
${OUTDIR}:
mkdir -p ${OUTDIR}
mkdir -p ${OUTDIR}
4 changes: 4 additions & 0 deletions papero.f90
Original file line number Diff line number Diff line change
Expand Up @@ -160,6 +160,10 @@ program papero
timeflag = .true.
#endif

#ifdef SPEC
specflag = .true.
#endif

call gcurv

errorcode = 1
Expand Down
13 changes: 11 additions & 2 deletions param.f90
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,9 @@ module param
real, allocatable, dimension(:,:,:) :: d_UsolidT_dxj
logical, allocatable, dimension(:,:,:) :: solid_mask

logical :: timeflag

logical :: timeflag = .false.
logical :: specflag = .false.


!==========================================================
!******* Grid indices**************************************
Expand Down Expand Up @@ -255,3 +256,11 @@ module mls_local
real, dimension(:,:,:), allocatable :: for_xc, for_yc, for_zc, for_temp
end module mls_local

module modspec
implicit none
integer*8 :: specplan
complex, allocatable :: uhat(:,:,:)
!complex, allocatable :: uhat(:)

end module modspec

6 changes: 3 additions & 3 deletions part.in
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
IMLSFOR IMLSSTR IMELT
1 1 0
0 0 0
rhop/rhof
10.0
1.0
GTSFX RAD
unitSphere_N4.gts 0.1
unitSphere_N3.gts 0.1
IREMESH PERCT_ETHRESH V_ON_VE_THRESH
0 0.7 27.0
3 changes: 3 additions & 0 deletions quit.f90
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,9 @@ subroutine QuitRoutine(tin,normalexit,errorcode)

call mem_dealloc
call dealloc_trigeo

if (specflag) call dealloc_spec

call FinalizeMPI

endif
Expand Down
169 changes: 169 additions & 0 deletions spectra.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,169 @@
! This subroutine is for computing 1D spectra on-the-fly during run-time

! For now, we only compute Exx(kx) (average taken in x and z), requiring only an fft in the x-direction
subroutine compute_1d_spectra
use local_arrays, only: vx
use modspec
use param
use mpih
use mpi_param, only: kstart, kend
implicit none
real, allocatable, dimension(:) :: Exx_kx
!real, allocatable, dimension(:,:,:) :: Exx

integer i,j,k
real :: dkx

allocate(Exx_kx(1:n1m/2+1))
!allocate(Exx(1:n1m/2+1,1:n2m, kstart:kend ) )

dkx = 2.0 * pi / xlen

Exx_kx = 0.0
! Exx = 0.0

do j = 1,n2m
do k = kstart, kend
call dfftw_execute_dft_r2c(specplan, vx(1:n1m,j,k), uhat(:,j,k))
!call dfftw_execute_dft(specplan, vx(:,j,k), uhat(:,j,k))

uhat(:,j,k) = uhat(:,j,k) / dble(n1m)
!Exx_kx(:) = Exx_kx(:) + ( real(uhat)**2 + aimag(uhat)**2 )
!Exx(:,j,k) = ( real(uhat(:,j,k))**2 + aimag(uhat(:,j,k))**2 )

Exx_kx = Exx_kx + ( real(uhat(:,j,k))**2 + aimag(uhat(:,j,k))**2 )

enddo
enddo

! do j = 1,n2m
! do k = kstart, kend
! Exx_kx(:) = Exx_kx(:) + Exx(:,j,k)
! enddo
! enddo

! do j = 1,n2m
! do k = kstart,kend
! !call dfftw_execute_dft(specplan, vx(:,j,k), uhat(:,j,k))
! !Exx_kx(i) = Exx_kx(i) + abs( uhat(i,j,k) )**2

! call dfftw_execute_dft(specplan, vx(:,j,k), uhat)
! uhat = uhat / dble(n1m)
! Exx_kx = Exx_kx + ( real(uhat)**2 + aimag(uhat)**2 )

! enddo
! enddo

call MPI_ALLREDUCE(MPI_IN_PLACE,Exx_kx,n1m/2+1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)

Exx_kx = Exx_kx / dble(n2m * n3m) / dkx

! Write to file
call writeSpec(Exx_kx , n1m/2+1 )

deallocate( Exx_kx )
!deallocate( Exx )

end subroutine compute_1d_spectra


! Initialisation of working memory, FFTW plans etc used in computing spectra
subroutine initSpectra
use local_arrays, only: vx
use modspec
use param
use mpi_param, only: kstart, kend
implicit none
integer :: my_n3m
integer FFTW_EXHAUSTIVE
parameter(FFTW_EXHAUSTIVE=64)

my_n3m = kend - kstart + 1
allocate(uhat(1:(n1m/2+1), 1:n2m, kstart:kend))
!allocate( uhat(1:(n1m/2+1) ) )


!call dfftw_plan_many_dft_r2c(specplan,1, n1m, n2m*my_n3m, vx(:,:,kstart:kend), [1, n2m*my_n3m], n2m*my_n3m, 1, uhat, [1, n2m*my_n3m], n2m*my_n3m, 1, FFTW_EXHAUSTIVE)
!call dfftw_plan_many_dft_r2c(specplan,1, n1m, n2m*my_n3m, vx(:,:,kstart:kend), n1m, 1, n1m, uhat, [1, n2m*my_n3m], n2m*my_n3m, 1, FFTW_EXHAUSTIVE)



call dfftw_plan_dft_r2c_1d(specplan,n1m,vx(1:n1m,1,kstart),uhat(:,1,kstart),FFTW_EXHAUSTIVE)

!call dfftw_plan_dft_r2c_1d(specplan,n1m,vx(:,1,kstart),uhat(:),FFTW_EXHAUSTIVE)


end subroutine initSpectra

subroutine dealloc_spec
use modspec
implicit none

call dfftw_destroy_plan(specplan)

if(allocated(uhat)) deallocate(uhat)

end subroutine dealloc_spec


subroutine writeSpec(Exx_kx, size_kx)
use param
use mpih
use mls_param
use hdf5

implicit none

character(70) filename
character(30) dataset
character(30) dsetname
character(5) ipfi
integer :: itime
real :: tprfi
integer :: rank = 1 ! Dataset rank
integer(hid_t) :: file_id ! File identifier
integer(hid_t) :: dset_id ! Dataset identifier
integer(hid_t) :: dspace_id ! Dataspace identifier
integer(hsize_t) :: dims(1)
integer :: error ! Error flag
logical :: fileexists


integer :: size_kx
real, dimension(size_kx) :: Exx_kx


tprfi = 1/tframe
itime=nint(time*tprfi)
write(ipfi,82)itime
82 format(i5.5)

filename ='spectra/Exx_kx_'//ipfi//'.h5'


if (myid.eq.0) then
call hdf5_create_blank_file(filename)
dsetname = trim("Exx_kx")
!call hdf_write_1d(Exx_kx,n1m/2+1,dataset)

!hdf_write_1d(dataset,d,dsetname)

dims(1) = n1m/2+1

call h5fopen_f(trim(filename), H5F_ACC_RDWR_F, file_id, error)
call h5screate_simple_f(rank, dims, dspace_id, error)

call h5lexists_f(file_id,dsetname,fileexists,error)
if(fileexists) call h5ldelete_f(file_id,dsetname,error)

call h5dcreate_f(file_id, dsetname, H5T_NATIVE_DOUBLE, dspace_id, dset_id, error)

call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, Exx_kx, dims, error)

call h5dclose_f(dset_id, error)
call h5sclose_f(dspace_id, error)
call h5fclose_f(file_id, error)

end if

end subroutine writeSpec

0 comments on commit d50199e

Please sign in to comment.