Skip to content

Commit

Permalink
fft
Browse files Browse the repository at this point in the history
  • Loading branch information
Anand committed Feb 19, 2025
1 parent e8c84d7 commit 8d79869
Show file tree
Hide file tree
Showing 4 changed files with 69 additions and 66 deletions.
4 changes: 4 additions & 0 deletions src/post_process/m_global_parameters.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,6 +217,8 @@ module m_global_parameters
integer :: flux_lim
logical, dimension(3) :: flux_wrt
logical :: E_wrt
logical :: fft_wrt
logical :: enstrophy_wrt
logical :: pres_wrt
logical, dimension(num_fluids_max) :: alpha_wrt
logical :: gamma_wrt
Expand Down Expand Up @@ -383,6 +385,8 @@ contains
parallel_io = .false.
file_per_process = .false.
E_wrt = .false.
fft_wrt = .false.
enstrophy_wrt = .false.
pres_wrt = .false.
alpha_wrt = .false.
gamma_wrt = .false.
Expand Down
9 changes: 8 additions & 1 deletion src/post_process/m_mpi_proxy.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,7 @@ contains
#:for VAR in [ 'cyl_coord', 'mpp_lim', 'mixture_err', &
& 'alt_soundspeed', 'hypoelasticity', 'parallel_io', 'rho_wrt', &
& 'E_wrt', 'pres_wrt', 'gamma_wrt', 'sim_data', &
& 'E_wrt', 'enstrophy_wrt', 'fft_wrt','pres_wrt', 'gamma_wrt', 'sim_data', &
& 'heat_ratio_wrt', 'pi_inf_wrt', 'pres_inf_wrt', 'cons_vars_wrt', &
& 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt', 'bubbles_euler', 'qbmm', &
& 'polytropic', 'polydisperse', 'file_per_process', 'relax', 'cf_wrt', &
Expand Down Expand Up @@ -359,6 +359,13 @@ contains
end if
!! Particular decomposition needed for FFT
if(fft_wrt) then
num_procs_z = num_procs
num_procs_x = 1
num_procs_y = 1
end if
! Checking whether the decomposition of the computational
! domain was successful
if (proc_rank == 0 .and. ierr == -1) then
Expand Down
120 changes: 55 additions & 65 deletions src/post_process/m_start_up.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ subroutine s_read_input_file
hypoelasticity, G, &
chem_wrt_Y, chem_wrt_T, avg_state, &
alpha_rho_wrt, rho_wrt, mom_wrt, vel_wrt, &
E_wrt, pres_wrt, alpha_wrt, gamma_wrt, &
E_wrt, fft_wrt, enstrophy_wrt, pres_wrt, alpha_wrt, gamma_wrt, &
heat_ratio_wrt, pi_inf_wrt, pres_inf_wrt, &
cons_vars_wrt, prim_vars_wrt, c_wrt, &
omega_wrt, qm_wrt, schlieren_wrt, schlieren_alpha, &
Expand Down Expand Up @@ -357,7 +357,17 @@ subroutine s_save_data(t_step, varname, pres, c, H)
! Adding the energy to the formatted database file
if (E_wrt .or. cons_vars_wrt) then

if(all((/bc_x%beg,bc_x%end,bc_y%beg,bc_y%end,bc_z%beg,bc_z%end/) == -1)) then
q_sf = q_cons_vf(E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
write (varname, '(A)') 'E'
call s_write_variable_to_formatted_database_file(varname, t_step)

varname(:) = ' '

end if

!Adding Energy cascade FFT
if(fft_wrt) then
if(num_procs == 1 .and. all((/bc_x%beg,bc_x%end,bc_y%beg,bc_y%end,bc_z%beg,bc_z%end/) == -1)) then

data_real = CMPLX(q_cons_vf(mom_idx%beg)%sf(0:m, 0:n, 0:p) / q_cons_vf(1)%sf(0:m, 0:n, 0:p), 0d0)

Expand Down Expand Up @@ -505,13 +515,6 @@ subroutine s_save_data(t_step, varname, pres, c, H)
end if
end do
end if

q_sf = q_cons_vf(E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
write (varname, '(A)') 'E'
call s_write_variable_to_formatted_database_file(varname, t_step)

varname(:) = ' '

end if

! Adding the elastic shear stresses to the formatted database file
Expand Down Expand Up @@ -660,57 +663,44 @@ subroutine s_save_data(t_step, varname, pres, c, H)
if (omega_wrt(i)) then

call s_derive_vorticity_component(i, q_prim_vf, q_sf)

if(i == 1) then

omega_tot = 0d0
do l = 0, p
do k = 0, n
do j = 0, m
omega_tot = omega_tot + q_cons_vf(1)%sf(j, k, l)*q_sf(j,k,l)**2d0
if(enstrophy_wrt) then
if(i == 1) then

omega_tot = 0d0
do l = 0, p
do k = 0, n
do j = 0, m
omega_tot = omega_tot + q_cons_vf(1)%sf(j, k, l)*q_sf(j,k,l)**2d0
end do
end do
end do
end do

if(proc_rank == 0) then
!print *, "OMEGA1", omega_tot, proc_rank
end if

else if(i == 2) then
do l = 0, p
do k = 0, n
do j = 0, m
omega_tot = omega_tot + q_cons_vf(1)%sf(j, k, l)*q_sf(j,k,l)**2d0
else if(i == 2) then
do l = 0, p
do k = 0, n
do j = 0, m
omega_tot = omega_tot + q_cons_vf(1)%sf(j, k, l)*q_sf(j,k,l)**2d0
end do
end do
end do
end do

if(proc_rank == 0) then
!print *, "OMEGA2", omega_tot, proc_rank
end if

else
do l = 0, p
do k = 0, n
do j = 0, m
omega_tot = omega_tot + q_cons_vf(1)%sf(j, k, l)*q_sf(j,k,l)**2d0
end do
else
do l = 0, p
do k = 0, n
do j = 0, m
omega_tot = omega_tot + q_cons_vf(1)%sf(j, k, l)*q_sf(j,k,l)**2d0
end do
end do
end do
end do

if(proc_rank == 0) then
!print *, "OMEGA3", omega_tot, proc_rank
end if
write(filename,'(a,i0,a)') 'omega_tot',proc_rank,'.dat'
inquire (FILE=filename, EXIST=file_exists)
if (file_exists) then
open (1, file=filename, position='append', status='old')
write (1, *) omega_tot, proc_rank, t_step
close (1)
else
open (1, file=filename, status='new')
write (1, *) omega_tot, proc_rank, t_step
close (1)
write(filename,'(a,i0,a)') 'omega_tot',proc_rank,'.dat'
inquire (FILE=filename, EXIST=file_exists)
if (file_exists) then
open (1, file=filename, position='append', status='old')
write (1, *) omega_tot, proc_rank, t_step
close (1)
else
open (1, file=filename, status='new')
write (1, *) omega_tot, proc_rank, t_step
close (1)
end if
end if
end if

Expand Down Expand Up @@ -870,18 +860,18 @@ subroutine s_initialize_modules
s_read_data_files => s_read_serial_data_files
else
s_read_data_files => s_read_parallel_data_files
end if
fftw_real_data = fftw_alloc_complex(int((m+1)*(n+1)*(p+1),c_size_t))
fftw_cmplx_data = fftw_alloc_complex(int((m+1)*(n+1)*(p+1),c_size_t))

call c_f_pointer(fftw_real_data, data_real, [m+1,n+1,p+1])
call c_f_pointer(fftw_cmplx_data, data_cmplx, [m+1,n+1,p+1])

fwd_plan = fftw_plan_dft_3d(m+1, n+1, p+1, data_real, data_cmplx, FFTW_FORWARD, FFTW_ESTIMATE)
end if

if(E_wrt) then
if(fft_wrt) then
allocate(vel1_real(1:m+1, 1:n+1, 1:p+1), vel2_real(1:m+1, 1:n+1, 1:p+1), vel3_real(1:m+1, 1:n+1, 1:p+1) , En_real(1:m+1, 1:n+1, 1:p+1))
allocate(En(1:m+1))
allocate(En(1:m+1))
if(num_procs == 1) then
fftw_real_data = fftw_alloc_complex(int((m+1)*(n+1)*(p+1),c_size_t))
fftw_cmplx_data = fftw_alloc_complex(int((m+1)*(n+1)*(p+1),c_size_t))
call c_f_pointer(fftw_real_data, data_real, [m+1,n+1,p+1])
call c_f_pointer(fftw_cmplx_data, data_cmplx, [m+1,n+1,p+1])
fwd_plan = fftw_plan_dft_3d(p+1, n+1, m+1, data_real, data_cmplx, FFTW_FORWARD, FFTW_ESTIMATE)
end if
end if
end subroutine s_initialize_modules

Expand Down
2 changes: 2 additions & 0 deletions toolchain/mfc/run/case_dicts.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -391,6 +391,8 @@ def analytic(self):
'flux_lim': ParamType.INT,
'flux_wrt': ParamType.LOG,
'E_wrt': ParamType.LOG,
'enstrophy_wrt': ParamType.LOG,
'fft_wrt': ParamType.LOG,
'pres_wrt': ParamType.LOG,
'alpha_wrt': ParamType.LOG,
'kappa_wrt': ParamType.LOG,
Expand Down

0 comments on commit 8d79869

Please sign in to comment.