Skip to content
Snippets Groups Projects
Commit e5cd9674 authored by Jean-Matthieu Etancelin's avatar Jean-Matthieu Etancelin
Browse files

Add patch to convert Parmes into SinglePrecision (to use: patch -p0 <... 

Add patch to convert Parmes into SinglePrecision (to use: patch -p0 < ParmesToSinglePrecision.patch)
parent bcc6a5b6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
HySoP/ParmesToSinglePrecision.patch 0 → 100644
+ 950
− 0
View file @ e5cd9674
diff --git parmepy/constants.py parmepy/constants.py
index c92cc42..8d85877 100644
--- parmepy/constants.py
+++ parmepy/constants.py
@@ -18,13 +18,13 @@ else:
PI = math.pi
# Set default type for real and integer numbers
-PARMES_REAL = np.float64
+PARMES_REAL = np.float32
# type for array indices
PARMES_INDEX = np.uint32
# type for integers
-PARMES_INTEGER = np.int64
+PARMES_INTEGER = np.int32
# float type for MPI messages
-PARMES_MPI_REAL = MPI.DOUBLE
+PARMES_MPI_REAL = MPI.REAL
## default array layout (fortran or C convention)
ORDER = 'F'
# to check array ordering with :
diff --git parmepy/f2py/fftw2py.f90 parmepy/f2py/fftw2py.f90
index 8645a9f..718eb19 100755
--- parmepy/f2py/fftw2py.f90
+++ parmepy/f2py/fftw2py.f90
@@ -5,7 +5,7 @@
module fftw2py
use client_data
- use parmesparam
+ use parmesparam_sp
!> 2d case
use fft2d
!> 3d case
@@ -81,7 +81,7 @@ contains
subroutine solve_poisson_2d(omega,velocity_x,velocity_y)
real(pk),dimension(:,:),intent(in):: omega
real(pk),dimension(size(omega,1),size(omega,2)),intent(out) :: velocity_x,velocity_y
- real(pk) :: start
+ real(8) :: start
!f2py intent(in,out) :: velocity_x,velocity_y
start = MPI_WTime()
@@ -117,7 +117,7 @@ contains
real(pk),dimension(:,:,:),intent(in):: omega_x,omega_y,omega_z
real(pk),dimension(size(omega_x,1),size(omega_y,2),size(omega_z,3)),intent(out) :: velocity_x,velocity_y,velocity_z
integer, dimension(3), intent(in) :: ghosts_vort, ghosts_velo
- real(pk) :: start
+ real(8) :: start
!f2py intent(in,out) :: velocity_x,velocity_y,velocity_z
start = MPI_WTime()
call r2c_3d(omega_x,omega_y,omega_z, ghosts_vort)
diff --git parmepy/f2py/scales2py.f90 parmepy/f2py/scales2py.f90
index e574369..391e127 100755
--- parmepy/f2py/scales2py.f90
+++ parmepy/f2py/scales2py.f90
@@ -5,7 +5,7 @@ use cart_topology, only : cart_create,set_group_size,discretisation_create,N_pro
use advec, only : advec_init,advec_step,advec_step_Inter_basic,advec_step_Inter_Two
use interpolation_velo, only : interpol_init
use mpi
-use parmesparam
+use parmesparam_sp
implicit none
@@ -91,7 +91,7 @@ contains
real(pk), dimension(size(vx,1),size(vx,2),size(vx,3)), intent(inout) :: scal
!f2py real(pk) intent(in,out), depend(size(vx,1)) :: scal
- real(pk) :: t0
+ real(8) :: t0
t0 = MPI_Wtime()
call advec_step(dt,vx,vy,vz,scal)
diff --git parmepy/operator/tests/test_velocity_correction.py parmepy/operator/tests/test_velocity_correction.py
old mode 100755
new mode 100644
diff --git setup.py.in setup.py.in
index 3c75c97..4a8c17a 100644
--- setup.py.in
+++ setup.py.in
@@ -71,8 +71,8 @@ if(enable_fortran is "ON"):
fortran_src.append(fortran_dir+'fftw2py.f90')
fftwdir = '@FFTWLIB@'
fftwdir = os.path.split(fftwdir)[0]
- parmeslib.append('fftw3')
- parmeslib.append('fftw3_mpi')
+ parmeslib.append('fftw3f')
+ parmeslib.append('fftw3f_mpi')
parmes_libdir.append(fftwdir)
withscales = '@WITH_SCALES@'
diff --git src/client_data.f90 src/client_data.f90
index 46b5268..77178d9 100755
--- src/client_data.f90
+++ src/client_data.f90
@@ -1,14 +1,14 @@
!> Some global parameters and variables
module client_data
- use MPI, only : MPI_DOUBLE_PRECISION
+ use MPI, only : MPI_REAL
use, intrinsic :: iso_c_binding ! required for fftw
implicit none
!> kind for real variables (simple or double precision)
- integer, parameter :: mk = kind(1.0d0) ! double precision
+ integer, parameter :: mk = kind(1.0) ! single precision
!> kind for real variables in mpi routines
- integer, parameter :: mpi_mk = MPI_DOUBLE_PRECISION
+ integer, parameter :: mpi_mk = MPI_REAL
!> Problem dimension (model, required for ppm to work properly)
integer, parameter :: dime = 2
!> Real dimension
@@ -26,7 +26,7 @@ module client_data
!> to activate (or not) screen output
logical,parameter :: verbose = .True.
!> i (sqrt(-1) ...)
- complex(C_DOUBLE_COMPLEX), parameter :: Icmplx = cmplx(0._mk,1._mk, kind=mk)
+ complex(C_FLOAT_COMPLEX), parameter :: Icmplx = cmplx(0._mk,1._mk, kind=mk)
!> tolerance used to compute error
real(mk), parameter :: tolerance = 1e-12
diff --git src/fftw/Poisson.f90 src/fftw/Poisson.f90
index 78b355c..8804e3d 100755
--- src/fftw/Poisson.f90
+++ src/fftw/Poisson.f90
@@ -53,7 +53,7 @@ contains
real(mk),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
real(mk),dimension(:,:,:),intent(inout) :: velocity_x,velocity_y,velocity_z
integer, dimension(3), intent(in) :: ghosts_w, ghosts_v
- real(mk) :: start
+ real(8) :: start
!! Compute fftw forward transform
!! Omega is used to initialize the fftw buffer for input field.
@@ -74,7 +74,7 @@ contains
real(mk),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
real(mk),dimension(:,:,:),intent(inout) :: velocity_x,velocity_y,velocity_z
- real(mk) :: start
+ real(8) :: start
!! Compute fftw forward transform
!! Omega is used to initialize the fftw buffer for input field.
diff --git src/fftw/fft2d.f90 src/fftw/fft2d.f90
index c56412f..7edbfd3 100755
--- src/fftw/fft2d.f90
+++ src/fftw/fft2d.f90
@@ -36,15 +36,15 @@ module fft2d
type(C_PTR) :: cbuffer2
!! Note Franck : check if local declarations of datain/out works and improve perfs.
!> Field (complex values) for fftw input
- complex(C_DOUBLE_COMPLEX), pointer :: datain1(:,:),datain2(:,:)
+ complex(C_FLOAT_COMPLEX), pointer :: datain1(:,:),datain2(:,:)
!> Field (real values) for fftw input
- real(C_DOUBLE), pointer :: rdatain1(:,:)
+ real(C_FLOAT), pointer :: rdatain1(:,:)
!> Field (complex values) for fftw (forward) output
- complex(C_DOUBLE_COMPLEX), pointer :: dataout1(:,:)
+ complex(C_FLOAT_COMPLEX), pointer :: dataout1(:,:)
!> Field (real values) for fftw output
- real(C_DOUBLE), pointer :: rdatain2(:,:)
+ real(C_FLOAT), pointer :: rdatain2(:,:)
!> Field (complex values) for fftw (forward) output
- complex(C_DOUBLE_COMPLEX), pointer :: dataout2(:,:)
+ complex(C_FLOAT_COMPLEX), pointer :: dataout2(:,:)
!> GLOBAL number of points in each direction
integer(C_INTPTR_T),pointer :: fft_resolution(:)
!> LOCAL resolution
@@ -52,13 +52,13 @@ module fft2d
!> Offset in the direction of distribution
integer(c_INTPTR_T),dimension(2) :: local_offset
!> wave numbers for fft in x direction
- real(C_DOUBLE), pointer :: kx(:)
+ real(C_FLOAT), pointer :: kx(:)
!> wave numbers for fft in y direction
- real(C_DOUBLE), pointer :: ky(:)
+ real(C_FLOAT), pointer :: ky(:)
!> log file for fftw
character(len=20),parameter :: filename ="parmesfftw.log"
!> normalization factor
- real(C_DOUBLE) :: normFFT
+ real(C_FLOAT) :: normFFT
!> true if all the allocation stuff for global variables has been done.
logical :: is2DUpToDate = .false.
@@ -81,7 +81,7 @@ contains
if(is2DUpToDate) return
! init fftw mpi context
- call fftw_mpi_init()
+ call fftwf_mpi_init()
if(rank==0) open(unit=21,file=filename,form="formatted")
@@ -90,27 +90,27 @@ contains
fft_resolution = resolution-1
! compute "optimal" size (according to fftw) for local date (warning : dimension reversal)
- alloc_local = fftw_mpi_local_size_2d_transposed(fft_resolution(c_Y),fft_resolution(c_X),main_comm,&
+ alloc_local = fftwf_mpi_local_size_2d_transposed(fft_resolution(c_Y),fft_resolution(c_X),main_comm,&
local_resolution(c_Y),local_offset(c_Y),local_resolution(c_X),local_offset(c_X));
! allocate local buffer (used to save datain/dataout1 ==> in-place transform!!)
- cbuffer1 = fftw_alloc_complex(alloc_local)
+ cbuffer1 = fftwf_alloc_complex(alloc_local)
! link datain and dataout1 to cbuffer, setting the right dimensions for each
call c_f_pointer(cbuffer1, datain1, [fft_resolution(c_X),local_resolution(c_Y)])
call c_f_pointer(cbuffer1, dataout1, [fft_resolution(c_Y),local_resolution(c_X)])
! second buffer used for backward transform. Used to copy dataout1 into dataout2 (input for backward transform and filter)
! and to save (in-place) the transform of the second component of the velocity
- cbuffer2 = fftw_alloc_complex(alloc_local)
+ cbuffer2 = fftwf_alloc_complex(alloc_local)
call c_f_pointer(cbuffer2, datain2, [fft_resolution(c_X),local_resolution(c_Y)])
call c_f_pointer(cbuffer2, dataout2, [fft_resolution(c_Y),local_resolution(c_X)])
! create MPI plan for in-place forward/backward DFT (note dimension reversal)
- plan_forward1 = fftw_mpi_plan_dft_2d(fft_resolution(c_Y), fft_resolution(c_X),datain1,dataout1,&
+ plan_forward1 = fftwf_mpi_plan_dft_2d(fft_resolution(c_Y), fft_resolution(c_X),datain1,dataout1,&
main_comm,FFTW_FORWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
- plan_backward1 = fftw_mpi_plan_dft_2d(fft_resolution(c_Y),fft_resolution(c_X),dataout1,datain1,&
+ plan_backward1 = fftwf_mpi_plan_dft_2d(fft_resolution(c_Y),fft_resolution(c_X),dataout1,datain1,&
main_comm,FFTW_BACKWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
- plan_backward2 = fftw_mpi_plan_dft_2d(fft_resolution(c_Y),fft_resolution(c_X),dataout2,datain2,&
+ plan_backward2 = fftwf_mpi_plan_dft_2d(fft_resolution(c_Y),fft_resolution(c_X),dataout2,datain2,&
main_comm,FFTW_BACKWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
call computeKxC(lengths(c_X))
@@ -139,7 +139,7 @@ contains
end do
! compute transform (as many times as desired)
- call fftw_mpi_execute_dft(plan_forward1, datain1, dataout1)
+ call fftwf_mpi_execute_dft(plan_forward1, datain1, dataout1)
!!$ do i = 1, fft_resolution(c_Y)
!!$ write(*,'(a,i5,a,16f10.4)') 'out[',rank,'] ', dataout1(i,1:local_resolution(c_X))
@@ -147,8 +147,8 @@ contains
!!$
call filter_poisson_2d()
- call fftw_mpi_execute_dft(plan_backward1, dataout1, datain1)
- call fftw_mpi_execute_dft(plan_backward2,dataout2,datain2)
+ call fftwf_mpi_execute_dft(plan_backward1, dataout1, datain1)
+ call fftwf_mpi_execute_dft(plan_backward2,dataout2,datain2)
do j = 1, local_resolution(c_Y)
do i = 1, fft_resolution(c_X)
velocity_x(i,j) = datain1(i,j)*normFFT
@@ -183,7 +183,7 @@ contains
if(is2DUpToDate) return
! init fftw mpi context
- call fftw_mpi_init()
+ call fftwf_mpi_init()
if(rank==0) open(unit=21,file=filename,form="formatted")
@@ -191,11 +191,11 @@ contains
fft_resolution(:) = resolution(:)-1
halfLength = fft_resolution(c_X)/2+1
! allocate local buffer (used to save datain/dataout1 ==> in-place transform!!)
- alloc_local = fftw_mpi_local_size_2d_transposed(fft_resolution(c_Y),halfLength,main_comm,local_resolution(c_Y),&
+ alloc_local = fftwf_mpi_local_size_2d_transposed(fft_resolution(c_Y),halfLength,main_comm,local_resolution(c_Y),&
local_offset(c_Y),local_resolution(c_X),local_offset(c_X));
! allocate local buffer (used to save datain/dataout1 ==> in-place transform!!)
- cbuffer1 = fftw_alloc_complex(alloc_local)
+ cbuffer1 = fftwf_alloc_complex(alloc_local)
! link rdatain1 and dataout1 to cbuffer, setting the right dimensions for each
call c_f_pointer(cbuffer1, rdatain1, [2*halfLength,local_resolution(c_Y)])
@@ -203,17 +203,17 @@ contains
! second buffer used for backward transform. Used to copy dataout1 into dataout2 (input for backward transform and filter)
! and to save (in-place) the transform of the second component of the velocity
- cbuffer2 = fftw_alloc_complex(alloc_local)
+ cbuffer2 = fftwf_alloc_complex(alloc_local)
call c_f_pointer(cbuffer2, rdatain2, [2*halfLength,local_resolution(c_Y)])
call c_f_pointer(cbuffer2, dataout2, [fft_resolution(c_Y),local_resolution(c_X)])
! create MPI plans for in-place forward/backward DFT (note dimension reversal)
- plan_forward1 = fftw_mpi_plan_dft_r2c_2d(fft_resolution(c_Y), fft_resolution(c_X), rdatain1, dataout1, &
+ plan_forward1 = fftwf_mpi_plan_dft_r2c_2d(fft_resolution(c_Y), fft_resolution(c_X), rdatain1, dataout1, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
- plan_backward1 = fftw_mpi_plan_dft_c2r_2d(fft_resolution(c_Y), fft_resolution(c_X), dataout1, rdatain1, &
+ plan_backward1 = fftwf_mpi_plan_dft_c2r_2d(fft_resolution(c_Y), fft_resolution(c_X), dataout1, rdatain1, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
- plan_backward2 = fftw_mpi_plan_dft_c2r_2d(fft_resolution(c_Y), fft_resolution(c_X), dataout2, rdatain2, &
+ plan_backward2 = fftwf_mpi_plan_dft_c2r_2d(fft_resolution(c_Y), fft_resolution(c_X), dataout2, rdatain2, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
call computeKx(lengths(c_X))
@@ -248,7 +248,7 @@ contains
!!$ end do
!!$
! compute transform (as many times as desired)
- call fftw_mpi_execute_dft_r2c(plan_forward1, rdatain1, dataout1)
+ call fftwf_mpi_execute_dft_r2c(plan_forward1, rdatain1, dataout1)
!!$ do i = 1, fft_resolution(c_Y)
!!$ write(*,'(a,i5,a,16f10.4)') 'aaaa[',rank,'] ', dataout1(i,1:local_resolution(c_X))
@@ -261,8 +261,8 @@ contains
real(mk),dimension(:,:),intent(inout) :: velocity_x,velocity_y
integer(C_INTPTR_T) :: i, j
- call fftw_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
- call fftw_mpi_execute_dft_c2r(plan_backward2,dataout2,rdatain2)
+ call fftwf_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
+ call fftwf_mpi_execute_dft_c2r(plan_backward2,dataout2,rdatain2)
do j = 1, local_resolution(c_Y)
do i = 1, fft_resolution(c_X)
velocity_x(i,j) = rdatain1(i,j)*normFFT
@@ -288,7 +288,7 @@ contains
real(mk),dimension(:,:),intent(inout) :: omega
integer(C_INTPTR_T) :: i, j
- call fftw_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
+ call fftwf_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
do j = 1, local_resolution(c_Y)
do i = 1, fft_resolution(c_X)
omega(i,j) = rdatain1(i,j)*normFFT
@@ -390,7 +390,7 @@ contains
subroutine filter_poisson_2d()
integer(C_INTPTR_T) :: i, j
- complex(C_DOUBLE_COMPLEX) :: coeff
+ complex(C_FLOAT_COMPLEX) :: coeff
if(local_offset(c_X)==0) then
if(local_offset(c_Y) == 0) then
dataout1(1,1) = 0.0
@@ -437,9 +437,9 @@ contains
subroutine filter_diffusion_2d(nudt)
- real(C_DOUBLE), intent(in) :: nudt
+ real(C_FLOAT), intent(in) :: nudt
integer(C_INTPTR_T) :: i, j
- complex(C_DOUBLE_COMPLEX) :: coeff
+ complex(C_FLOAT_COMPLEX) :: coeff
do i = 1,local_resolution(c_X)
do j = 1, fft_resolution(c_Y)
@@ -452,20 +452,20 @@ contains
!> Clean fftw context (free memory, plans ...)
subroutine cleanFFTW_2d()
- call fftw_destroy_plan(plan_forward1)
- call fftw_destroy_plan(plan_backward1)
- !call fftw_destroy_plan(plan_forward2)
- !call fftw_destroy_plan(plan_backward2)
- call fftw_free(cbuffer1)
- call fftw_free(cbuffer2)
- call fftw_mpi_cleanup()
+ call fftwf_destroy_plan(plan_forward1)
+ call fftwf_destroy_plan(plan_backward1)
+ !call fftwf_destroy_plan(plan_forward2)
+ !call fftwf_destroy_plan(plan_backward2)
+ call fftwf_free(cbuffer1)
+ call fftwf_free(cbuffer2)
+ call fftwf_mpi_cleanup()
deallocate(fft_resolution)
if(rank==0) close(21)
end subroutine cleanFFTW_2d
subroutine fft2d_diagnostics(nbelem)
integer(C_INTPTR_T), intent(in) :: nbelem
- complex(C_DOUBLE_COMPLEX) :: memoryAllocated
+ complex(C_FLOAT_COMPLEX) :: memoryAllocated
memoryAllocated = real(nbelem*sizeof(memoryAllocated),mk)*1e-6
write(*,'(a,i5,a,i12,f10.2)') '[',rank,'] size of each buffer (elements / memory in MB):', &
nbelem, memoryAllocated
@@ -499,10 +499,10 @@ contains
integer(C_INTPTR_T), dimension(2) :: n
!> Field (real values) for fftw input
- real(C_DOUBLE), pointer :: rdatain1Many(:,:,:)
+ real(C_FLOAT), pointer :: rdatain1Many(:,:,:)
! init fftw mpi context
- call fftw_mpi_init()
+ call fftwf_mpi_init()
howmany = 1
if(rank==0) open(unit=21,file=filename,form="formatted")
@@ -514,12 +514,12 @@ contains
n(1) = fft_resolution(2)
n(2) = halfLength
! allocate local buffer (used to save datain/dataout1 ==> in-place transform!!)
- alloc_local = fftw_mpi_local_size_many_transposed(2,n,howmany,FFTW_MPI_DEFAULT_BLOCK,&
+ alloc_local = fftwf_mpi_local_size_many_transposed(2,n,howmany,FFTW_MPI_DEFAULT_BLOCK,&
FFTW_MPI_DEFAULT_BLOCK,main_comm,local_resolution(c_Y),&
local_offset(c_Y),local_resolution(c_X),local_offset(c_X));
! allocate local buffer (used to save datain/dataout1 ==> in-place transform!!)
- cbuffer1 = fftw_alloc_complex(alloc_local)
+ cbuffer1 = fftwf_alloc_complex(alloc_local)
! link rdatain1 and dataout1 to cbuffer, setting the right dimensions for each
call c_f_pointer(cbuffer1, rdatain1Many, [howmany,2*halfLength,local_resolution(c_Y)])
@@ -527,17 +527,17 @@ contains
! second buffer used for backward transform. Used to copy dataout1 into dataout2 (input for backward transform and filter)
! and to save (in-place) the transform of the second component of the velocity
- cbuffer2 = fftw_alloc_complex(alloc_local)
+ cbuffer2 = fftwf_alloc_complex(alloc_local)
call c_f_pointer(cbuffer2, rdatain1Many, [howmany,2*halfLength,local_resolution(c_Y)])
call c_f_pointer(cbuffer2, dataout2, [fft_resolution(c_Y),local_resolution(c_X)])
! create MPI plans for in-place forward/backward DFT (note dimension reversal)
- plan_forward1 = fftw_mpi_plan_dft_r2c_2d(fft_resolution(c_Y), fft_resolution(c_X), rdatain1Many, dataout1, &
+ plan_forward1 = fftwf_mpi_plan_dft_r2c_2d(fft_resolution(c_Y), fft_resolution(c_X), rdatain1Many, dataout1, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
- plan_backward1 = fftw_mpi_plan_dft_c2r_2d(fft_resolution(c_Y), fft_resolution(c_X), dataout1, rdatain1, &
+ plan_backward1 = fftwf_mpi_plan_dft_c2r_2d(fft_resolution(c_Y), fft_resolution(c_X), dataout1, rdatain1, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
- plan_backward2 = fftw_mpi_plan_dft_c2r_2d(fft_resolution(c_Y), fft_resolution(c_X), dataout2, rdatain2, &
+ plan_backward2 = fftwf_mpi_plan_dft_c2r_2d(fft_resolution(c_Y), fft_resolution(c_X), dataout2, rdatain2, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
call computeKx(lengths(c_X))
diff --git src/fftw/fft3d.f90 src/fftw/fft3d.f90
index 1b06457..a857855 100755
--- src/fftw/fft3d.f90
+++ src/fftw/fft3d.f90
@@ -40,15 +40,15 @@ module fft3d
type(C_PTR) :: cbuffer3
!! Note Franck : check if local declarations of datain/out works and improve perfs.
!> Field (complex values) for fftw input
- complex(C_DOUBLE_COMPLEX), pointer :: datain1(:,:,:)=>NULL(), datain2(:,:,:)=>NULL(), datain3(:,:,:)=>NULL()
+ complex(C_FLOAT_COMPLEX), pointer :: datain1(:,:,:)=>NULL(), datain2(:,:,:)=>NULL(), datain3(:,:,:)=>NULL()
!> Field (real values) for fftw input (these are only pointers to the cbuffers)
- real(C_DOUBLE), pointer :: rdatain1(:,:,:)=>NULL() ,rdatain2(:,:,:)=>NULL() ,rdatain3(:,:,:)=>NULL()
+ real(C_FLOAT), pointer :: rdatain1(:,:,:)=>NULL() ,rdatain2(:,:,:)=>NULL() ,rdatain3(:,:,:)=>NULL()
!> Field (real values) for fftw input in the fftw-many case
- real(C_DOUBLE), pointer :: rdatain_many(:,:,:,:)=>NULL()
+ real(C_FLOAT), pointer :: rdatain_many(:,:,:,:)=>NULL()
!> Field (complex values) for fftw (forward) output
- complex(C_DOUBLE_COMPLEX), pointer :: dataout1(:,:,:)=>NULL() ,dataout2(:,:,:)=>NULL() ,dataout3(:,:,:)=>NULL()
+ complex(C_FLOAT_COMPLEX), pointer :: dataout1(:,:,:)=>NULL() ,dataout2(:,:,:)=>NULL() ,dataout3(:,:,:)=>NULL()
!> Field (complex values) for fftw (forward) output in the fftw-many case
- complex(C_DOUBLE_COMPLEX), pointer :: dataout_many(:,:,:,:)=>NULL()
+ complex(C_FLOAT_COMPLEX), pointer :: dataout_many(:,:,:,:)=>NULL()
!> GLOBAL number of points in each direction on which fft is applied (--> corresponds to "real" resolution - 1)
integer(C_INTPTR_T),pointer :: fft_resolution(:)=>NULL()
!> LOCAL number of points for fft
@@ -56,15 +56,15 @@ module fft3d
!> Offset in the direction of distribution
integer(c_INTPTR_T),dimension(3) :: local_offset
!> wave numbers for fft in x direction
- real(C_DOUBLE), pointer :: kx(:)
+ real(C_FLOAT), pointer :: kx(:)
!> wave numbers for fft in y direction
- real(C_DOUBLE), pointer :: ky(:)
+ real(C_FLOAT), pointer :: ky(:)
!> wave numbers for fft in z direction
- real(C_DOUBLE), pointer :: kz(:)
+ real(C_FLOAT), pointer :: kz(:)
!> log file for fftw
character(len=20),parameter :: filename ="parmesfftw.log"
!> normalization factor
- real(C_DOUBLE) :: normFFT
+ real(C_FLOAT) :: normFFT
!> true if we use fftw-many routines
logical :: manycase
!> true if all the allocation stuff for global variables has been done.
@@ -89,7 +89,7 @@ contains
if(is3DUpToDate) return
! init fftw mpi context
- call fftw_mpi_init()
+ call fftwf_mpi_init()
if(rank==0) open(unit=21,file=filename,form="formatted")
@@ -98,7 +98,7 @@ contains
fft_resolution(:) = resolution(:)-1
! compute "optimal" size (according to fftw) for local data (warning : dimension reversal)
- alloc_local = fftw_mpi_local_size_3d_transposed(fft_resolution(c_Z),fft_resolution(c_Y),fft_resolution(c_X),main_comm,&
+ alloc_local = fftwf_mpi_local_size_3d_transposed(fft_resolution(c_Z),fft_resolution(c_Y),fft_resolution(c_X),main_comm,&
local_resolution(c_Z),local_offset(c_Z),local_resolution(c_Y),local_offset(c_Y));
! Set a default value for c_X components.
@@ -106,35 +106,35 @@ contains
local_resolution(c_X) = fft_resolution(c_X)
! allocate local buffer (used to save datain/dataout ==> in-place transform!!)
- cbuffer1 = fftw_alloc_complex(alloc_local)
+ cbuffer1 = fftwf_alloc_complex(alloc_local)
! link datain and dataout to cbuffer, setting the right dimensions for each
call c_f_pointer(cbuffer1, datain1, [fft_resolution(c_X),fft_resolution(c_Y),local_resolution(c_Z)])
call c_f_pointer(cbuffer1, dataout1, [fft_resolution(c_X),fft_resolution(c_Z),local_resolution(c_Y)])
! second buffer used for backward transform. Used to copy dataout into dataout2 (input for backward transform and filter)
! and to save (in-place) the transform of the second component of the velocity
- cbuffer2 = fftw_alloc_complex(alloc_local)
+ cbuffer2 = fftwf_alloc_complex(alloc_local)
call c_f_pointer(cbuffer2, datain2, [fft_resolution(c_X),fft_resolution(c_Y),local_resolution(c_Z)])
call c_f_pointer(cbuffer2, dataout2, [fft_resolution(c_X),fft_resolution(c_Z),local_resolution(c_Y)])
! second buffer used for backward transform. Used to copy dataout into dataout2 (input for backward transform and filter)
! and to save (in-place) the transform of the second component of the velocity
- cbuffer3 = fftw_alloc_complex(alloc_local)
+ cbuffer3 = fftwf_alloc_complex(alloc_local)
call c_f_pointer(cbuffer3, datain3, [fft_resolution(c_X),fft_resolution(c_Y),local_resolution(c_Z)])
call c_f_pointer(cbuffer3, dataout3, [fft_resolution(c_X),fft_resolution(c_Z),local_resolution(c_Y)])
! create MPI plan for in-place forward/backward DFT (note dimension reversal)
- plan_forward1 = fftw_mpi_plan_dft_3d(fft_resolution(c_Z), fft_resolution(c_Y), fft_resolution(c_X),datain1,dataout1,&
+ plan_forward1 = fftwf_mpi_plan_dft_3d(fft_resolution(c_Z), fft_resolution(c_Y), fft_resolution(c_X),datain1,dataout1,&
main_comm,FFTW_FORWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
- plan_backward1 = fftw_mpi_plan_dft_3d(fft_resolution(c_Z),fft_resolution(c_Y),fft_resolution(c_X),dataout1,datain1,&
+ plan_backward1 = fftwf_mpi_plan_dft_3d(fft_resolution(c_Z),fft_resolution(c_Y),fft_resolution(c_X),dataout1,datain1,&
main_comm,FFTW_BACKWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
- plan_forward2 = fftw_mpi_plan_dft_3d(fft_resolution(c_Z), fft_resolution(c_Y), fft_resolution(c_X),datain2,dataout2,&
+ plan_forward2 = fftwf_mpi_plan_dft_3d(fft_resolution(c_Z), fft_resolution(c_Y), fft_resolution(c_X),datain2,dataout2,&
main_comm,FFTW_FORWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
- plan_backward2 = fftw_mpi_plan_dft_3d(fft_resolution(c_Z),fft_resolution(c_Y),fft_resolution(c_X),dataout2,datain2,&
+ plan_backward2 = fftwf_mpi_plan_dft_3d(fft_resolution(c_Z),fft_resolution(c_Y),fft_resolution(c_X),dataout2,datain2,&
main_comm,FFTW_BACKWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
- plan_forward3 = fftw_mpi_plan_dft_3d(fft_resolution(c_Z), fft_resolution(c_Y), fft_resolution(c_X),datain3,dataout3,&
+ plan_forward3 = fftwf_mpi_plan_dft_3d(fft_resolution(c_Z), fft_resolution(c_Y), fft_resolution(c_X),datain3,dataout3,&
main_comm,FFTW_FORWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
- plan_backward3 = fftw_mpi_plan_dft_3d(fft_resolution(c_Z),fft_resolution(c_Y),fft_resolution(c_X),dataout3,datain3,&
+ plan_backward3 = fftwf_mpi_plan_dft_3d(fft_resolution(c_Z),fft_resolution(c_Y),fft_resolution(c_X),dataout3,datain3,&
main_comm,FFTW_BACKWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
call computeKx(lengths(c_X))
@@ -173,17 +173,17 @@ contains
end do
end do
! compute transform (as many times as desired)
- call fftw_mpi_execute_dft(plan_forward1, datain1, dataout1)
- call fftw_mpi_execute_dft(plan_forward2, datain2, dataout2)
- call fftw_mpi_execute_dft(plan_forward3, datain3, dataout3)
+ call fftwf_mpi_execute_dft(plan_forward1, datain1, dataout1)
+ call fftwf_mpi_execute_dft(plan_forward2, datain2, dataout2)
+ call fftwf_mpi_execute_dft(plan_forward3, datain3, dataout3)
! apply poisson filter
call filter_poisson_3d()
! inverse transform to retrieve velocity
- call fftw_mpi_execute_dft(plan_backward1, dataout1,datain1)
- call fftw_mpi_execute_dft(plan_backward2,dataout2,datain2)
- call fftw_mpi_execute_dft(plan_backward3,dataout3,datain3)
+ call fftwf_mpi_execute_dft(plan_backward1, dataout1,datain1)
+ call fftwf_mpi_execute_dft(plan_backward2,dataout2,datain2)
+ call fftwf_mpi_execute_dft(plan_backward3,dataout3,datain3)
do k =1, local_resolution(c_Z)
do j = 1, fft_resolution(c_Y)
do i = 1, fft_resolution(c_X)
@@ -213,7 +213,7 @@ contains
if(is3DUpToDate) return
! init fftw mpi context
- call fftw_mpi_init()
+ call fftwf_mpi_init()
if(rank==0) open(unit=21,file=filename,form="formatted")
allocate(fft_resolution(3))
@@ -221,7 +221,7 @@ contains
halfLength = fft_resolution(c_X)/2+1
! compute "optimal" size (according to fftw) for local data (warning : dimension reversal)
- alloc_local = fftw_mpi_local_size_3d_transposed(fft_resolution(c_Z),fft_resolution(c_Y),halfLength,&
+ alloc_local = fftwf_mpi_local_size_3d_transposed(fft_resolution(c_Z),fft_resolution(c_Y),halfLength,&
main_comm,local_resolution(c_Z),local_offset(c_Z),local_resolution(c_Y),local_offset(c_Y));
! init c_X part. This is required to compute kx with the same function in 2d and 3d cases.
@@ -229,9 +229,9 @@ contains
local_resolution(c_X) = halfLength
! allocate local buffer (used to save datain/dataout ==> in-place transform!!)
- cbuffer1 = fftw_alloc_complex(alloc_local)
- cbuffer2 = fftw_alloc_complex(alloc_local)
- cbuffer3 = fftw_alloc_complex(alloc_local)
+ cbuffer1 = fftwf_alloc_complex(alloc_local)
+ cbuffer2 = fftwf_alloc_complex(alloc_local)
+ cbuffer3 = fftwf_alloc_complex(alloc_local)
! link rdatain and dataout to cbuffer, setting the right dimensions for each
call c_f_pointer(cbuffer1, rdatain1, [2*halfLength,fft_resolution(c_Y),local_resolution(c_Z)])
@@ -246,17 +246,17 @@ contains
rdatain3 = 0.0
! create MPI plans for in-place forward/backward DFT (note dimension reversal)
- plan_forward1 = fftw_mpi_plan_dft_r2c_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), rdatain1, dataout1, &
+ plan_forward1 = fftwf_mpi_plan_dft_r2c_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), rdatain1, dataout1, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
- plan_backward1 = fftw_mpi_plan_dft_c2r_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), dataout1, rdatain1, &
+ plan_backward1 = fftwf_mpi_plan_dft_c2r_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), dataout1, rdatain1, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
- plan_forward2 = fftw_mpi_plan_dft_r2c_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), rdatain2, dataout2, &
+ plan_forward2 = fftwf_mpi_plan_dft_r2c_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), rdatain2, dataout2, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
- plan_backward2 = fftw_mpi_plan_dft_c2r_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), dataout2, rdatain2, &
+ plan_backward2 = fftwf_mpi_plan_dft_c2r_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), dataout2, rdatain2, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
- plan_forward3 = fftw_mpi_plan_dft_r2c_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), rdatain3, dataout3, &
+ plan_forward3 = fftwf_mpi_plan_dft_r2c_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), rdatain3, dataout3, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
- plan_backward3 = fftw_mpi_plan_dft_c2r_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), dataout3, rdatain3, &
+ plan_backward3 = fftwf_mpi_plan_dft_c2r_3d(fft_resolution(c_Z),fft_resolution(c_Y), fft_resolution(c_X), dataout3, rdatain3, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
call computeKx(lengths(c_X))
@@ -279,7 +279,7 @@ contains
real(mk),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
integer, dimension(3), intent(in) :: ghosts
- real(mk) :: start
+ real(8) :: start
integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
! ig, jg, kg are used to take into account
@@ -301,9 +301,9 @@ contains
! compute transforms for each component
start = MPI_WTIME()
- call fftw_mpi_execute_dft_r2c(plan_forward1, rdatain1, dataout1)
- call fftw_mpi_execute_dft_r2c(plan_forward2, rdatain2, dataout2)
- call fftw_mpi_execute_dft_r2c(plan_forward3, rdatain3, dataout3)
+ call fftwf_mpi_execute_dft_r2c(plan_forward1, rdatain1, dataout1)
+ call fftwf_mpi_execute_dft_r2c(plan_forward2, rdatain2, dataout2)
+ call fftwf_mpi_execute_dft_r2c(plan_forward3, rdatain3, dataout3)
!!print *, "r2c time = ", MPI_WTIME() - start
end subroutine r2c_3d
@@ -316,12 +316,12 @@ contains
subroutine c2r_3d(velocity_x,velocity_y,velocity_z, ghosts)
real(mk),dimension(:,:,:),intent(inout) :: velocity_x,velocity_y,velocity_z
integer, dimension(3), intent(in) :: ghosts
- real(mk) :: start
+ real(8) :: start
integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
start = MPI_WTIME()
- call fftw_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
- call fftw_mpi_execute_dft_c2r(plan_backward2,dataout2,rdatain2)
- call fftw_mpi_execute_dft_c2r(plan_backward3,dataout3,rdatain3)
+ call fftwf_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
+ call fftwf_mpi_execute_dft_c2r(plan_backward2,dataout2,rdatain2)
+ call fftwf_mpi_execute_dft_c2r(plan_backward3,dataout3,rdatain3)
!! print *, "c2r time : ", MPI_WTIME() -start
! copy back to velocity and normalisation
do k =1, local_resolution(c_Z)
@@ -346,7 +346,7 @@ contains
real(mk),dimension(:,:,:),intent(in) :: omega
integer, dimension(3), intent(in) :: ghosts
- real(mk) :: start
+ real(8) :: start
integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
! ig, jg, kg are used to take into account
@@ -366,7 +366,7 @@ contains
! compute transforms for each component
start = MPI_WTIME()
- call fftw_mpi_execute_dft_r2c(plan_forward1, rdatain1, dataout1)
+ call fftwf_mpi_execute_dft_r2c(plan_forward1, rdatain1, dataout1)
!!print *, "r2c time = ", MPI_WTIME() - start
end subroutine r2c_scalar_3d
@@ -377,10 +377,10 @@ contains
subroutine c2r_scalar_3d(omega, ghosts)
real(mk),dimension(:,:,:),intent(inout) :: omega
integer, dimension(3), intent(in) :: ghosts
- real(mk) :: start
+ real(8) :: start
integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
start = MPI_WTIME()
- call fftw_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
+ call fftwf_mpi_execute_dft_c2r(plan_backward1,dataout1,rdatain1)
!! print *, "c2r time : ", MPI_WTIME() -start
! copy back to velocity and normalisation
do k =1, local_resolution(c_Z)
@@ -412,7 +412,7 @@ contains
integer(C_INTPTR_T),dimension(3) :: n
! init fftw mpi context
- call fftw_mpi_init()
+ call fftwf_mpi_init()
blocksize = FFTW_MPI_DEFAULT_BLOCK
if(rank==0) open(unit=21,file=filename,form="formatted")
allocate(fft_resolution(3))
@@ -423,7 +423,7 @@ contains
n(3) = halfLength
howmany = 3
! compute "optimal" size (according to fftw) for local data (warning : dimension reversal)
- alloc_local = fftw_mpi_local_size_many_transposed(3,n,howmany,blocksize,blocksize,&
+ alloc_local = fftwf_mpi_local_size_many_transposed(3,n,howmany,blocksize,blocksize,&
main_comm,local_resolution(c_Z),local_offset(c_Z),local_resolution(c_Y),local_offset(c_Y));
! init c_X part. This is required to compute kx with the same function in 2d and 3d cases.
@@ -431,7 +431,7 @@ contains
local_resolution(c_X) = halfLength
! allocate local buffer (used to save datain/dataout ==> in-place transform!!)
- cbuffer1 = fftw_alloc_complex(alloc_local)
+ cbuffer1 = fftwf_alloc_complex(alloc_local)
! link rdatain and dataout to cbuffer, setting the right dimensions for each
call c_f_pointer(cbuffer1, rdatain_many, [howmany,2*halfLength,fft_resolution(c_Y),local_resolution(c_Z)])
@@ -440,9 +440,9 @@ contains
! create MPI plans for in-place forward/backward DFT (note dimension reversal)
n(3) = fft_resolution(c_X)
- plan_forward1 = fftw_mpi_plan_many_dft_r2c(3,n,howmany,blocksize,blocksize, rdatain_many, dataout_many, &
+ plan_forward1 = fftwf_mpi_plan_many_dft_r2c(3,n,howmany,blocksize,blocksize, rdatain_many, dataout_many, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
- plan_backward1 = fftw_mpi_plan_many_dft_c2r(3,n,howmany,blocksize,blocksize, dataout_many, rdatain_many, &
+ plan_backward1 = fftwf_mpi_plan_many_dft_c2r(3,n,howmany,blocksize,blocksize, dataout_many, rdatain_many, &
main_comm,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
call computeKx(lengths(c_X))
call computeKy(lengths(c_Y))
@@ -464,7 +464,7 @@ contains
subroutine r2c_3d_many(omega_x,omega_y,omega_z)
real(mk),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
- real(mk) :: start
+ real(8) :: start
integer(C_INTPTR_T) :: i,j,k
! init
@@ -480,7 +480,7 @@ contains
! compute transform (as many times as desired)
start = MPI_WTIME()
- call fftw_mpi_execute_dft_r2c(plan_forward1, rdatain_many, dataout_many)
+ call fftwf_mpi_execute_dft_r2c(plan_forward1, rdatain_many, dataout_many)
!! print *, "r2c time = ", MPI_WTIME() - start
end subroutine r2c_3d_many
@@ -491,11 +491,11 @@ contains
!! @param[in,out] velocity_z 3d scalar field, z-component of the output vector field
subroutine c2r_3d_many(velocity_x,velocity_y,velocity_z)
real(mk),dimension(:,:,:),intent(inout) :: velocity_x,velocity_y,velocity_z
- real(mk) :: start
+ real(8) :: start
integer(C_INTPTR_T) :: i,j,k
start = MPI_WTIME()
- call fftw_mpi_execute_dft_c2r(plan_backward1,dataout_many,rdatain_many)
+ call fftwf_mpi_execute_dft_c2r(plan_backward1,dataout_many,rdatain_many)
!! print *, "c2r time : ", MPI_WTIME() -start
do k =1, local_resolution(c_Z)
do j = 1, fft_resolution(c_Y)
@@ -536,7 +536,7 @@ contains
!> Computation of frequencies coeff, over distributed direction(s)
!> @param lengths size of the domain
subroutine computeKy(length)
- real(C_DOUBLE), intent(in) :: length
+ real(C_FLOAT), intent(in) :: length
!! Local loops indices
integer(C_INTPTR_T) :: i
@@ -589,8 +589,8 @@ contains
subroutine filter_poisson_3d()
integer(C_INTPTR_T) :: i,j,k
- complex(C_DOUBLE_COMPLEX) :: coeff
- complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2
+ complex(C_FLOAT_COMPLEX) :: coeff
+ complex(C_FLOAT_COMPLEX) :: buffer1,buffer2
! Set first coeff (check for "all freq = 0" case)
if(local_offset(c_Y) == 0) then
@@ -648,10 +648,10 @@ contains
!! @param[in] nudt \f$ \nu\times dt\f$, diffusion coefficient times current time step
subroutine filter_curl_diffusion_3d(nudt)
- real(C_DOUBLE), intent(in) :: nudt
+ real(C_FLOAT), intent(in) :: nudt
integer(C_INTPTR_T) :: i,j,k
- complex(C_DOUBLE_COMPLEX) :: coeff
- complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2
+ complex(C_FLOAT_COMPLEX) :: coeff
+ complex(C_FLOAT_COMPLEX) :: buffer1,buffer2
!! mind the transpose -> index inversion between y and z
do j = 1,local_resolution(c_Y)
@@ -674,9 +674,9 @@ contains
!! @param[in] nudt \f$ \nu\times dt\f$, diffusion coefficient times current time step
subroutine filter_diffusion_3d(nudt)
- real(C_DOUBLE), intent(in) :: nudt
+ real(C_FLOAT), intent(in) :: nudt
integer(C_INTPTR_T) :: i,j,k
- complex(C_DOUBLE_COMPLEX) :: coeff
+ complex(C_FLOAT_COMPLEX) :: coeff
!! mind the transpose -> index inversion between y and z
do j = 1,local_resolution(c_Y)
@@ -697,8 +697,8 @@ contains
subroutine filter_curl_3d()
integer(C_INTPTR_T) :: i,j,k
- complex(C_DOUBLE_COMPLEX) :: coeff
- complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2
+ complex(C_FLOAT_COMPLEX) :: coeff
+ complex(C_FLOAT_COMPLEX) :: buffer1,buffer2
!! mind the transpose -> index inversion between y and z
do j = 1,local_resolution(c_Y)
@@ -721,8 +721,8 @@ contains
subroutine filter_projection_om_3d()
integer(C_INTPTR_T) :: i,j,k
- complex(C_DOUBLE_COMPLEX) :: coeff
- complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2,buffer3
+ complex(C_FLOAT_COMPLEX) :: coeff
+ complex(C_FLOAT_COMPLEX) :: buffer1,buffer2,buffer3
! Set first coeff (check for "all freq = 0" case)
if(local_offset(c_Y) == 0) then
@@ -783,9 +783,9 @@ contains
!! @param[in] dxf, dyf, dzf: grid filter size = domainLength/(CoarseRes-1)
subroutine filter_multires_om_3d(dxf, dyf, dzf)
- real(C_DOUBLE), intent(in) :: dxf, dyf, dzf
+ real(C_FLOAT), intent(in) :: dxf, dyf, dzf
integer(C_INTPTR_T) :: i,j,k
- real(C_DOUBLE) :: kxc, kyc, kzc
+ real(C_FLOAT) :: kxc, kyc, kzc
kxc = pi / dxf
kyc = pi / dyf
@@ -810,7 +810,7 @@ contains
!! pressure from velocity in the Fourier space
subroutine filter_pressure_3d()
integer(C_INTPTR_T) :: i,j,k
- complex(C_DOUBLE_COMPLEX) :: coeff
+ complex(C_FLOAT_COMPLEX) :: coeff
! Set first coeff (check for "all freq = 0" case)
if(local_offset(c_Y) == 0) then
@@ -851,8 +851,8 @@ contains
subroutine filter_poisson_3d_many()
integer(C_INTPTR_T) :: i,j,k
- complex(C_DOUBLE_COMPLEX) :: coeff
- complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2
+ complex(C_FLOAT_COMPLEX) :: coeff
+ complex(C_FLOAT_COMPLEX) :: buffer1,buffer2
! Set first coeff (check for "all freq = 0" case)
if(local_offset(c_Y) == 0) then
@@ -908,10 +908,10 @@ contains
!! @param[in] nudt \f$ \nu\times dt\f$, diffusion coefficient times current time step
subroutine filter_diffusion_3d_many(nudt)
- real(C_DOUBLE), intent(in) :: nudt
+ real(C_FLOAT), intent(in) :: nudt
integer(C_INTPTR_T) :: i,j,k
- complex(C_DOUBLE_COMPLEX) :: coeff
- complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2
+ complex(C_FLOAT_COMPLEX) :: coeff
+ complex(C_FLOAT_COMPLEX) :: buffer1,buffer2
!! mind the transpose -> index inversion between y and z
do j = 1,local_resolution(c_Y)
@@ -931,18 +931,18 @@ contains
!> Clean fftw context (free memory, plans ...)
subroutine cleanFFTW_3d()
- call fftw_destroy_plan(plan_forward1)
- call fftw_destroy_plan(plan_backward1)
+ call fftwf_destroy_plan(plan_forward1)
+ call fftwf_destroy_plan(plan_backward1)
if(.not.manycase) then
- call fftw_destroy_plan(plan_forward2)
- call fftw_destroy_plan(plan_backward2)
- call fftw_destroy_plan(plan_forward3)
- call fftw_destroy_plan(plan_backward3)
- call fftw_free(cbuffer2)
- call fftw_free(cbuffer3)
+ call fftwf_destroy_plan(plan_forward2)
+ call fftwf_destroy_plan(plan_backward2)
+ call fftwf_destroy_plan(plan_forward3)
+ call fftwf_destroy_plan(plan_backward3)
+ call fftwf_free(cbuffer2)
+ call fftwf_free(cbuffer3)
endif
- call fftw_free(cbuffer1)
- call fftw_mpi_cleanup()
+ call fftwf_free(cbuffer1)
+ call fftwf_mpi_cleanup()
deallocate(fft_resolution)
deallocate(kx,ky,kz)
if(rank==0) close(21)
@@ -953,7 +953,7 @@ contains
integer(C_INTPTR_T), intent(in) :: nbelem
! number of buffers used for fftw
integer, optional,intent(in) :: howmany
- complex(C_DOUBLE_COMPLEX) :: memoryAllocated
+ complex(C_FLOAT_COMPLEX) :: memoryAllocated
integer :: nbFields
if(present(howmany)) then
diff --git src/main/main.f90 src/main/main.f90
index 38cc6a5..0fc8f5f 100755
--- src/main/main.f90
+++ src/main/main.f90
@@ -9,7 +9,7 @@ use vectorcalculus
implicit none
integer :: info
-real(mk) :: start, end
+real(8) :: start, end
!complex(mk), dimension(resolution(1),resolution(2)) :: omega,velocity_x,velocity_y
call MPI_Init(info)
@@ -116,7 +116,8 @@ contains
real(mk),dimension(3) :: lengths,step
integer, dimension(3) :: ghosts_v, ghosts_w
integer(C_INTPTR_T),dimension(3) :: nfft,offset
- real(mk) :: error,start
+ real(mk) :: error
+ real(8) :: start
logical :: ok
if (rank==0) print *, " ======= Test 3D Poisson (r2c) solver for resolution ", resolution
diff --git src/scalesInterface/precision_tools.f90 src/scalesInterface/precision_tools.f90
index 5dacf3f..4ae1a6a 100644
--- src/scalesInterface/precision_tools.f90
+++ src/scalesInterface/precision_tools.f90
@@ -20,15 +20,15 @@
!------------------------------------------------------------------------------
MODULE precision_tools
- use mpi, only: MPI_DOUBLE_PRECISION
+ use mpi, only: MPI_REAL
implicit None
!> Floats precision
INTEGER, PARAMETER :: SP = kind(1.0)
INTEGER, PARAMETER :: DP = kind(1.0d0)
- INTEGER, PARAMETER :: WP = DP
+ INTEGER, PARAMETER :: WP = SP
!> the MPI type for REAL exchanges in simple or double precision
- INTEGER, parameter :: MPI_REAL_WP = MPI_DOUBLE_PRECISION
+ INTEGER, parameter :: MPI_REAL_WP = MPI_REAL
REAL(WP), PRIVATE :: sample_real_at_WP
REAL(WP), PARAMETER :: MAX_REAL_WP = HUGE(sample_real_at_WP)
INTEGER, PRIVATE :: sample_int
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment