From ff3a4e10500748aa74be0810b42e58f84611634a Mon Sep 17 00:00:00 2001
From: Jean-Matthieu Etancelin <jean-matthieu.etancelin@imag.fr>
Date: Fri, 31 Jan 2014 13:43:15 +0000
Subject: [PATCH] Fix real precision in SCALES (removing implicit conversions).
Fix bug in computing remeshing indices from particle adimentional position
(integer translation after floor position).
---
HySoP/hysop/f2py/fftw2py.f90 | 90 ++++----
HySoP/hysop/f2py/parameters.f90 | 11 +
HySoP/hysop/f2py/scales2py.f90 | 5 +-
HySoP/hysop/operator/energy_enstrophy.py | 4 +-
HySoP/src/VectorCalc.f90 | 60 +++---
HySoP/src/client_data.f90 | 16 +-
HySoP/src/fftw/fft2d.f90 | 120 +++++------
HySoP/src/fftw/fft3d.f90 | 204 +++++++++---------
HySoP/src/scalesInterface/particles/advec.f90 | 4 +-
.../src/scalesInterface/particles/advecX.f90 | 16 +-
.../src/scalesInterface/particles/advecY.f90 | 16 +-
.../src/scalesInterface/particles/advecZ.f90 | 16 +-
.../particles/advec_common_interpol.F90 | 16 +-
.../particles/advec_common_remesh.F90 | 8 +-
.../advec_line/advec_common_line.f90 | 78 +++----
.../advec_line/advec_remesh_formula_line.f90 | 96 ++++-----
.../advec_line/advec_remesh_line.F90 | 16 +-
.../particles/advec_remesh_Mprime.f90 | 90 +++++---
.../particles/advec_remesh_lambda.f90 | 144 ++++++-------
.../particles/interpolation_velo.f90 | 8 +-
HySoP/src/scalesInterface/precision_tools.f90 | 2 +-
21 files changed, 540 insertions(+), 480 deletions(-)
diff --git a/HySoP/hysop/f2py/fftw2py.f90 b/HySoP/hysop/f2py/fftw2py.f90
index edf1c4a99..6e290e6e5 100755
--- a/HySoP/hysop/f2py/fftw2py.f90
+++ b/HySoP/hysop/f2py/fftw2py.f90
@@ -1,4 +1,4 @@
-!> @file fftw2py.f90
+!> @file fftw2py.f90
!! Fortran to python interface file.
!> Interface to mpi-fftw (fortran) utilities
@@ -13,7 +13,7 @@ module fftw2py
implicit none
contains
-
+
!> Initialisation of fftw context : create plans and memory buffers
!! @param[in] resolution global resolution of the discrete domain
!! @param[in] lengths width of each side of the domain
@@ -21,7 +21,7 @@ contains
!! @param[out] offset absolute index of the first component of the local field
subroutine init_fftw_solver(resolution,lengths,datashape,offset,dim,fftw_type_real)
- integer, intent(in) :: dim
+ integer, intent(in) :: dim
integer, dimension(dim),intent(in) :: resolution
real(pk),dimension(dim), intent(in) :: lengths
integer(ik), dimension(dim), intent(out) :: datashape
@@ -29,10 +29,10 @@ contains
logical, optional :: fftw_type_real
!f2py optional :: dim=len(resolution)
!f2py intent(hide) dim
- !f2py logical optional, intent(in) :: fftw_type_real = 1
-
+ !f2py logical optional, intent(in) :: fftw_type_real = 1
+
if(fftw_type_real) then
- if(dim == 2) then
+ if(dim == 2) then
print*, "Init fftw/poisson solver for a 2d problem"
call init_r2c_2d(resolution,lengths)
else
@@ -40,7 +40,7 @@ contains
call init_r2c_3d(resolution,lengths)
end if
else
- if(dim == 2) then
+ if(dim == 2) then
print*, "Init fftw/poisson solver for a 2d problem"
call init_c2c_2d(resolution,lengths)
else
@@ -59,17 +59,17 @@ contains
!> Free memory allocated for fftw-related objects (plans and buffers)
subroutine clean_fftw_solver(dim)
- integer, intent(in) :: dim
- if(dim == 2) then
+ integer, intent(in) :: dim
+ if(dim == 2) then
call cleanFFTW_2d()
else
call cleanFFTW_3d()
end if
end subroutine clean_fftw_solver
-
- !> Solve
+
+ !> Solve
!! \f[ \nabla (\nabla \times velocity) = - \omega \f]
- !! velocity being a 2D vector field and omega a 2D scalar field.
+ !! velocity being a 2D vector field and omega a 2D scalar field.
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
@@ -78,31 +78,31 @@ contains
start = MPI_WTime()
call r2c_2d(omega)
-
+
call filter_poisson_2d()
call c2r_2d(velocity_x,velocity_y)
!!print *, "fortran resolution time : ", MPI_WTime() - start
-
+
end subroutine solve_poisson_2d
- !> Solve
+ !> Solve
!! \f{eqnarray*} \frac{\partial \omega}{\partial t} &=& \nu \Delta \omega \f}
- !! omega being a 2D scalar field.
+ !! omega being a 2D scalar field.
subroutine solve_diffusion_2d(nudt, omega)
real(pk), intent(in) :: nudt
real(pk),dimension(:,:),intent(inout):: omega
!f2py intent(in,out) :: omega
call r2c_2d(omega)
-
+
call filter_diffusion_2d(nudt)
call c2r_scalar_2d(omega)
-
+
end subroutine solve_diffusion_2d
- !> Solve
+ !> Solve
!! \f{eqnarray*} \Delta \psi &=& - \omega \\ velocity = \nabla\times\psi \f}
!! velocity and omega being 3D vector fields.
subroutine solve_poisson_3d(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_z, ghosts_vort, ghosts_velo)
@@ -110,30 +110,30 @@ contains
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
- !f2py intent(in,out) :: velocity_x,velocity_y,velocity_z
+ !f2py intent(in,out) :: velocity_x,velocity_y,velocity_z
start = MPI_WTime()
call r2c_3d(omega_x,omega_y,omega_z, ghosts_vort)
-
+
call filter_poisson_3d()
-
+
call c2r_3d(velocity_x,velocity_y,velocity_z, ghosts_velo)
!!print *, "fortran resolution time : ", MPI_WTime() - start
end subroutine solve_poisson_3d
- !> Solve
+ !> Solve
!! \f{eqnarray*} \Delta \psi &=& - \omega \\ velocity = \nabla\times\psi \f}
- !! velocity being a 2D complex vector field and omega a 2D complex scalar field.
+ !! velocity being a 2D complex vector field and omega a 2D complex scalar field.
subroutine solve_poisson_2d_c(omega,velocity_x,velocity_y)
complex(pk),dimension(:,:),intent(in):: omega
complex(pk),dimension(size(omega,1),size(omega,2)),intent(out) :: velocity_x,velocity_y
!f2py intent(in,out) :: velocity_x,velocity_y
call c2c_2d(omega,velocity_x,velocity_y)
-
+
end subroutine solve_poisson_2d_c
- !> Solve
+ !> Solve
!! \f{eqnarray*} \Delta \psi &=& - \omega \\ velocity = \nabla\times\psi \f}
!! velocity and omega being 3D complex vector fields.
subroutine solve_poisson_3d_c(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_Z)
@@ -142,10 +142,10 @@ contains
!f2py intent(in,out) :: velocity_x,velocity_y,velocity_z
call c2c_3d(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_z)
-
+
end subroutine solve_poisson_3d_c
- !> Solve
+ !> Solve
!! \f{eqnarray*} \omega &=& \nabla \times v \\ \frac{\partial \omega}{\partial t} &=& \nu \Delta \omega \f}
!! velocity and omega being 3D vector fields.
subroutine solve_curl_diffusion_3d(nudt,velocity_x,velocity_y,velocity_z,omega_x,omega_y,omega_z, ghosts_velo, ghosts_vort)
@@ -154,16 +154,16 @@ contains
real(pk),dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(out) :: omega_x,omega_y,omega_z
integer, dimension(3), intent(in) :: ghosts_vort, ghosts_velo
!f2py intent(in,out) :: omega_x,omega_y,omega_z
-
+
call r2c_3d(velocity_x,velocity_y,velocity_z, ghosts_velo)
-
+
call filter_curl_diffusion_3d(nudt)
-
+
call c2r_3d(omega_x,omega_y,omega_z, ghosts_vort)
end subroutine solve_curl_diffusion_3d
- !> Solve
+ !> Solve
!! \f{eqnarray*} \frac{\partial \omega}{\partial t} &=& \nu \Delta \omega \f}
!! omega being 3D vector field.
subroutine solve_diffusion_3d(nudt,omega_x,omega_y,omega_z, ghosts)
@@ -171,11 +171,11 @@ contains
real(pk),dimension(:,:,:),intent(inout):: omega_x,omega_y,omega_z
integer, dimension(3), intent(in) :: ghosts
!f2py intent(in,out) :: omega_x,omega_y,omega_z
-
+
call r2c_3d(omega_x,omega_y,omega_z, ghosts)
-
+
call filter_diffusion_3d(nudt)
-
+
call c2r_3d(omega_x,omega_y,omega_z, ghosts)
end subroutine solve_diffusion_3d
@@ -189,9 +189,9 @@ contains
!f2py intent(in,out) :: omega_x,omega_y,omega_z
call r2c_3d(omega_x,omega_y,omega_z, ghosts)
-
+
call filter_projection_om_3d()
-
+
call c2r_3d(omega_x,omega_y,omega_z, ghosts)
end subroutine projection_om_3d
@@ -207,9 +207,9 @@ contains
!f2py intent(in,out) :: omega_x,omega_y,omega_z
call r2c_3d(omega_x,omega_y,omega_z, ghosts)
-
+
call filter_multires_om_3d(dxf, dyf, dzf)
-
+
call c2r_3d(omega_x,omega_y,omega_z, ghosts)
end subroutine multires_om_3d
@@ -229,14 +229,14 @@ contains
!f2py intent(in,out) :: omega_x,omega_y,omega_z
call r2c_scalar_3d(rhs, ghosts)
-
+
call filter_pressure_3d()
-
+
call c2r_scalar_3d(pressure, ghosts)
end subroutine pressure_3d
- !> Solve
+ !> Solve
!! \f{eqnarray*} \omega &=& \nabla \times v
!! velocity and omega being 3D vector fields.
subroutine solve_curl_3d(velocity_x,velocity_y,velocity_z,omega_x,omega_y,omega_z, ghosts_velo, ghosts_vort)
@@ -244,11 +244,11 @@ contains
real(pk),dimension(size(velocity_x,1),size(velocity_y,2),size(velocity_z,3)),intent(out) :: omega_x,omega_y,omega_z
integer, dimension(3), intent(in) :: ghosts_velo, ghosts_vort
!f2py intent(in,out) :: omega_x,omega_y,omega_z
-
+
call r2c_3d(velocity_x,velocity_y,velocity_z, ghosts_velo)
-
+
call filter_curl_3d()
-
+
call c2r_3d(omega_x,omega_y,omega_z, ghosts_vort)
end subroutine solve_curl_3d
diff --git a/HySoP/hysop/f2py/parameters.f90 b/HySoP/hysop/f2py/parameters.f90
index 60f1a9165..bb281b27a 100755
--- a/HySoP/hysop/f2py/parameters.f90
+++ b/HySoP/hysop/f2py/parameters.f90
@@ -12,3 +12,14 @@ module parmesparam
integer, parameter :: ik = 8
end module parmesparam
+
+module parmesparam_sp
+
+ implicit none
+
+ ! double precision kind
+ integer, parameter :: pk = 4
+ ! integer precision kind
+ integer, parameter :: ik = 4
+
+end module parmesparam_sp
diff --git a/HySoP/hysop/f2py/scales2py.f90 b/HySoP/hysop/f2py/scales2py.f90
index e392eb384..62b469c8b 100755
--- a/HySoP/hysop/f2py/scales2py.f90
+++ b/HySoP/hysop/f2py/scales2py.f90
@@ -7,6 +7,7 @@ use interpolation_velo, only : interpol_init
use mpi
use parmesparam
+
implicit none
contains
@@ -26,7 +27,7 @@ contains
integer(ik), dimension(dim), intent(out) :: offset
character(len=*), optional, intent(in) :: order, dim_split
real(pk), optional, intent(out) :: stab_coeff
- !f2py optional , depend(ncells) :: dim=len(ncells)
+ !f2py integer optional , depend(ncells) :: dim=len(ncells)
!f2py intent(hide) dim
!f2py character(*) optional, intent(in) :: order = 'p_O2'
!f2py, depends(dim), intent(in) :: topodims
@@ -86,7 +87,7 @@ contains
real(pk), intent(in) :: dt
real(pk), dimension(:,:,:), intent(in) :: vx, vy, vz
real(pk), dimension(size(vx,1),size(vx,2),size(vx,3)), intent(inout) :: scal
- !f2py intent(in,out), depend(size(vx,1)) :: scal
+ !f2py real(pk) intent(in,out), depend(size(vx,1)) :: scal
real(pk) :: t0
diff --git a/HySoP/hysop/operator/energy_enstrophy.py b/HySoP/hysop/operator/energy_enstrophy.py
index 6126c8d14..3a66362c1 100644
--- a/HySoP/hysop/operator/energy_enstrophy.py
+++ b/HySoP/hysop/operator/energy_enstrophy.py
@@ -29,7 +29,7 @@ class Energy_enstrophy(Monitoring):
@param velocity field
@param vorticity field
@param viscosity : kinematic viscosity
- @param isNormalized : boolean indicating whether the enstrophy
+ @param isNormalized : boolean indicating whether the enstrophy
and energy values have to be normalized by the domain lengths.
@param topo : the topology on which we want to monitor the fields
@param frequency : output file producing frequency.
@@ -188,6 +188,6 @@ class Energy_enstrophy(Monitoring):
## nu_effEnstrophy))
def finalize(self):
- pass
+ pass
# if self._topo.rank == 0:
# self.f.close()
diff --git a/HySoP/src/VectorCalc.f90 b/HySoP/src/VectorCalc.f90
index a3117af24..098d9b5a3 100755
--- a/HySoP/src/VectorCalc.f90
+++ b/HySoP/src/VectorCalc.f90
@@ -1,6 +1,6 @@
!> routines to compute curl, cross prod ...
!! WARNING : many of the following routines are out-of-date with
-!! remaining "ppm" type call of field shapes.
+!! remaining "ppm" type call of field shapes.
!! \todo : clean everything and move to python.
module vectorcalculus
@@ -13,27 +13,27 @@ module vectorcalculus
integer, parameter,private :: nsublist = 1
contains
-
+
!> compute \f[ fieldout = \nabla \times fieldin \f] using 4th order finite differences
subroutine curlDF4(fieldin,fieldout,resolution,step)
!> input field
real(mk), dimension(:,:,:,:,:), pointer :: fieldin
- !> output field
+ !> output field
real(mk), dimension(:,:,:,:,:), pointer :: fieldout
!> the local resolution
integer, dimension(dim3),intent(in) :: resolution
!> size of mesh step in each dir
real(mk), dimension(dim3),intent(in) :: step
-
+
real(mk) :: facx, facy, facz
integer :: i,j,k
fieldout = 0.0
-
+
facx = 1.0/(12.0*step(c_X))
facy = 1.0/(12.0*step(c_Y))
facz = 1.0/(12.0*step(c_Z))
-
+
do k=1,resolution(c_Z)
do j=1,resolution(c_Y)
do i=1,resolution(c_X)
@@ -59,7 +59,7 @@ contains
enddo
enddo
enddo
-
+
end subroutine curlDF4
!> Computes strech and diffusion terms. This is a copy of Adrien's code.
@@ -152,7 +152,7 @@ contains
!diffusion
!----------
- tx(c_X)= - &
+ tx(c_X)= - &
vorticity(c_X,i+2,j,k,nsublist) + 16.*&
vorticity(c_X,i+1,j,k,nsublist) - 30.*&
vorticity(c_X,i,j,k,nsublist) + 16.*&
@@ -311,7 +311,7 @@ contains
end do
end subroutine computeStretch
-
+
function iscloseto_3d(x,y,tol,step,error)
logical :: iscloseto_3d
@@ -322,15 +322,15 @@ contains
real(mk) :: res
integer :: info
real(mk) :: h3
- res = sum(abs(x-y)**2)
+ res = sum(abs(x-y)**2)
iscloseto_3d = .false.
call MPI_AllReduce(res,error,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info)
h3 = product(step(:))
error = sqrt(h3*error)
if(error < tol) iscloseto_3d = .True.
-
+
end function iscloseto_3d
-
+
function iscloseto_2d(x,y,tol,step,error)
logical :: iscloseto_2d
@@ -341,13 +341,13 @@ contains
real(mk) :: res
integer :: info
real(mk) :: h3
- res = sum(abs(x-y)**2)
+ res = sum(abs(x-y)**2)
iscloseto_2d = .false.
call MPI_AllReduce(res,error,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info)
h3 = product(step(:))
error = sqrt(h3*error)
if(error < tol) iscloseto_2d = .True.
-
+
end function iscloseto_2d
function iscloseto_3dc(x,y,tol,step,error)
@@ -359,16 +359,16 @@ contains
real(mk) :: res
integer :: info
real(mk) :: h3
- res = sum(abs(real(x)-real(y))**2)
- res = res + sum(abs(aimag(x)-aimag(y))**2)
+ res = sum(abs(real(x,mk)-real(y,mk))**2)
+ res = res + sum(abs(aimag(x)-aimag(y))**2)
iscloseto_3dc = .false.
call MPI_AllReduce(res,error,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info)
h3 = product(step(:))
error = sqrt(h3*error)
if(error < tol) iscloseto_3dc = .True.
-
+
end function iscloseto_3dc
-
+
function iscloseto_2dc(x,y,tol,step,error)
logical :: iscloseto_2dc
@@ -379,19 +379,19 @@ contains
real(mk) :: res
integer :: info
real(mk) :: h3
- res = sum(abs(real(x)-real(y))**2)
- res = res + sum(abs(aimag(x)-aimag(y))**2)
+ res = sum(abs(real(x,mk)-real(y,mk))**2)
+ res = res + sum(abs(aimag(x)-aimag(y))**2)
iscloseto_2dc = .false.
call MPI_AllReduce(res,error,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,info)
h3 = product(step(:))
error = sqrt(h3*error)
if(error < tol) iscloseto_2dc = .True.
-
+
end function iscloseto_2dc
-
+
!!$ !> compute euclidian norm of a 3D real field
!!$ function norm2_f3d(field,resolution,step)
-!!$
+!!$
!!$ real(mk), dimension(3) :: norm2
!!$ real(mk), dimension(:,:,:), pointer :: field
!!$ !> the local resolution
@@ -412,7 +412,7 @@ contains
!!$ call MPI_Reduce(buffer,norm2,dime,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,info)
!!$ if(rank == 0) norm2 = sqrt(h3*norm2)
!!$ end function norm2
-!!$
+!!$
!!$ function norm22d(field,resolution,step)
!!$ real(mk), dimension(2) :: norm22d
!!$ real(mk), dimension(:,:), pointer :: field
@@ -433,9 +433,9 @@ contains
!!$ ! Norm is computed only on proc 0
!!$ call MPI_Reduce(buffer,norm22d,dime,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,info)
!!$ if(rank == 0) norm22d = sqrt(h2*norm22d)
-!!$
+!!$
!!$ end function norm22d
-!!$
+!!$
!!$ function norm22dC(field,resolution,step)
!!$ real(mk), dimension(2) :: norm22d
!!$ complex(mk), dimension(:,:), pointer :: field
@@ -456,7 +456,7 @@ contains
!!$ ! Norm is computed only on proc 0
!!$ call MPI_Reduce(buffer,norm22d,dime,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,info)
!!$ if(rank == 0) norm22d = sqrt(h2*norm22d)
-!!$
+!!$
!!$ end function norm22dC
!!$
!!$ function normInf(field,resolution)
@@ -469,12 +469,12 @@ contains
!!$ integer :: i,info
!!$ do i = 1,dime
!!$ buffer(i) = maxval(abs(field(i,1:resolution(c_X),1:resolution(c_Y),1:resolution(c_Z))))
-!!$ end do
+!!$ end do
!!$ ! Norm is computed only on proc 0
!!$ call MPI_Reduce(buffer,normInf,dime,MPI_DOUBLE_PRECISION,MPI_MAX,0,MPI_COMM_WORLD,info)
!!$
!!$ end function normInf
-
+
function cross_prod(v1,v2)
real(mk), dimension(3), intent(in) :: v1
real(mk), dimension(3), intent(in) :: v2
@@ -497,7 +497,7 @@ contains
!> mesh step size
real(mk),dimension(3),intent(in) :: step
- ! --- nabla vect_component ---
+ ! --- nabla vect_component ---
nabla(c_X) = (vect(i+1,j,k) - vect(i-1,j,k))/(2.*step(c_X)) ! d/dx u_component
nabla(c_Y) = (vect(i,j+1,k) - vect(i,j-1,k))/(2.*step(c_Y)) ! d/dy u_component
nabla(c_Z) = (vect(i,j,k+1) - vect(i,j,k-1))/(2.*step(c_Z)) ! d/dz u_component
diff --git a/HySoP/src/client_data.f90 b/HySoP/src/client_data.f90
index 1c6944ad3..ba2280bfb 100755
--- a/HySoP/src/client_data.f90
+++ b/HySoP/src/client_data.f90
@@ -3,11 +3,11 @@ module client_data
use MPI, only : MPI_DOUBLE_PRECISION, MPI_REAL,MPI_COMM_WORLD
use, intrinsic :: iso_c_binding ! required for fftw
- implicit none
-
+ implicit none
+
!> kind for real variables (simple or double precision)
- integer, parameter :: mk = kind(1.0d0) ! double precision
- !integer, parameter :: mk = kind(1.0) ! simple precision
+ integer, parameter :: mk = kind(1.0d0) ! double precision
+ !integer, parameter :: mk = kind(1.0) ! simple precision
!> kind for real variables in mpi routines
integer, parameter :: mpi_mk = MPI_DOUBLE_PRECISION
!integer, parameter :: mpi_mk = MPI_REAL
@@ -17,11 +17,11 @@ module client_data
integer, parameter :: dim3 = 3
!> Pi constant
real(mk), parameter :: pi = 4.0*atan(1.0_mk)
- !> Rank of the mpi current process
+ !> Rank of the mpi current process
integer :: rank ! current mpi-processus rank
!> Total number of mpi process
- integer :: nbprocs
- !> trick to identify coordinates in a more user-friendly way
+ integer :: nbprocs
+ !> trick to identify coordinates in a more user-friendly way
integer,parameter :: c_X=1,c_Y=2,c_Z=3
!> to activate (or not) screen output
logical,parameter :: verbose = .True.
@@ -29,5 +29,5 @@ module client_data
complex(C_DOUBLE_COMPLEX), parameter :: Icmplx = cmplx(0._mk,1._mk, kind=mk)
!> tolerance used to compute error
real(mk), parameter :: tolerance = 1e-12
-
+
end module client_data
diff --git a/HySoP/src/fftw/fft2d.f90 b/HySoP/src/fftw/fft2d.f90
index 7e5610f0e..d4a5bf26e 100755
--- a/HySoP/src/fftw/fft2d.f90
+++ b/HySoP/src/fftw/fft2d.f90
@@ -7,12 +7,12 @@
!! input fields which are real. The names of these routines contain "r2c".
!! \li 3 - fftw routines for the "real to complex" case : solves the problem for
!! input fields which are real and using the "many" interface of the fftw.
-!! It means that transforms are applied to the 2 input fields at the same time.
-!! Names of these routines contain "many".
-!!
-!! Obviously, all the above cases should lead to the same results. By default
+!! It means that transforms are applied to the 2 input fields at the same time.
+!! Names of these routines contain "many".
+!!
+!! Obviously, all the above cases should lead to the same results. By default
!! case 2 must be chosen (if input is real). Case 1 and 3 are more or less
-!! dedicated to tests and validation.
+!! dedicated to tests and validation.
module fft2d
use client_data
@@ -24,8 +24,8 @@ module fft2d
public :: init_c2c_2d,init_r2c_2d,c2c_2d,c2r_2d,c2r_scalar_2d,r2c_2d,cleanFFTW_2d,&
filter_poisson_2d,getParamatersTopologyFFTW2d,filter_diffusion_2d, init_r2c_2dBIS
-
-
+
+
!> plan for fftw "c2c" forward or r2c transform
type(C_PTR) :: plan_forward1, plan_forward2
!> plan for fftw "c2c" backward or c2r transform
@@ -39,15 +39,15 @@ module fft2d
complex(C_DOUBLE_COMPLEX), pointer :: datain1(:,:),datain2(:,:)
!> Field (real values) for fftw input
real(C_DOUBLE), pointer :: rdatain1(:,:)
- !> Field (complex values) for fftw (forward) output
+ !> Field (complex values) for fftw (forward) output
complex(C_DOUBLE_COMPLEX), pointer :: dataout1(:,:)
!> Field (real values) for fftw output
real(C_DOUBLE), pointer :: rdatain2(:,:)
- !> Field (complex values) for fftw (forward) output
+ !> Field (complex values) for fftw (forward) output
complex(C_DOUBLE_COMPLEX), pointer :: dataout2(:,:)
!> GLOBAL number of points in each direction
integer(C_INTPTR_T),pointer :: fft_resolution(:)
- !> LOCAL resolution
+ !> LOCAL resolution
integer(c_INTPTR_T),dimension(2) :: local_resolution
!> Offset in the direction of distribution
integer(c_INTPTR_T),dimension(2) :: local_offset
@@ -64,7 +64,7 @@ module fft2d
contains
!========================================================================
- ! Complex to complex transforms
+ ! Complex to complex transforms
!========================================================================
!> Initialisation of the fftw context for complex to complex transforms (forward and backward)
@@ -77,18 +77,18 @@ contains
!! Size of the local memory required for fftw (cbuffer)
integer(C_INTPTR_T) :: alloc_local
-
+
if(is2DUpToDate) return
! init fftw mpi context
call fftw_mpi_init()
-
+
if(rank==0) open(unit=21,file=filename,form="formatted")
-
+
! set fft resolution
allocate(fft_resolution(2))
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),MPI_COMM_WORLD,&
local_resolution(c_Y),local_offset(c_Y),local_resolution(c_X),local_offset(c_X));
@@ -112,7 +112,7 @@ contains
MPI_COMM_WORLD,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,&
MPI_COMM_WORLD,FFTW_BACKWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
-
+
call computeKxC(lengths(c_X))
call computeKy(lengths(c_Y))
normFFT = 1./(fft_resolution(c_X)*fft_resolution(c_Y))
@@ -129,9 +129,9 @@ contains
subroutine c2c_2d(inputData,velocity_x,velocity_y)
complex(mk),dimension(:,:) :: velocity_x,velocity_y
complex(mk),dimension(:,:),intent(in) :: inputData
-
+
integer(C_INTPTR_T) :: i, j
-
+
do j = 1, local_resolution(c_Y)
do i = 1, fft_resolution(c_X)
datain1(i, j) = inputData(i,j)
@@ -144,9 +144,9 @@ contains
!!$ do i = 1, fft_resolution(c_Y)
!!$ write(*,'(a,i5,a,16f10.4)') 'out[',rank,'] ', dataout1(i,1:local_resolution(c_X))
!!$ end do
-!!$
+!!$
call filter_poisson_2d()
-
+
call fftw_mpi_execute_dft(plan_backward1, dataout1, datain1)
call fftw_mpi_execute_dft(plan_backward2,dataout2,datain2)
do j = 1, local_resolution(c_Y)
@@ -155,7 +155,7 @@ contains
velocity_y(i,j) = datain2(i,j)*normFFT
end do
end do
-
+
!!$ do i = 1, fft_resolution(c_X)
!!$ write(*,'(a,i5,a,16f10.4)') 'vxx[',rank,'] ', velocity_x(i,1:local_resolution(c_Y))
!!$ end do
@@ -168,9 +168,9 @@ contains
end subroutine c2c_2d
!========================================================================
- ! Real to complex transforms
+ ! Real to complex transforms
!========================================================================
-
+
!> Initialisation of the fftw context for real to complex transforms (forward and backward)
!! @param[in] resolution global domain resolution
subroutine init_r2c_2d(resolution,lengths)
@@ -179,9 +179,9 @@ contains
real(mk),dimension(2), intent(in) :: lengths
!! Size of the local memory required for fftw (cbuffer)
integer(C_INTPTR_T) :: alloc_local,halfLength
-
+
if(is2DUpToDate) return
-
+
! init fftw mpi context
call fftw_mpi_init()
@@ -193,14 +193,14 @@ contains
! allocate local buffer (used to save datain/dataout1 ==> in-place transform!!)
alloc_local = fftw_mpi_local_size_2d_transposed(fft_resolution(c_Y),halfLength,MPI_COMM_WORLD,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)
! link rdatain1 and dataout1 to cbuffer, setting the right dimensions for each
call c_f_pointer(cbuffer1, rdatain1, [2*halfLength,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)
@@ -220,7 +220,7 @@ contains
call computeKy(lengths(c_Y))
normFFT = 1./(fft_resolution(c_X)*fft_resolution(c_Y))
call fft2d_diagnostics(alloc_local)
-!!$
+!!$
!!$ write(*,'(a,i5,a,16f10.4)') 'kx[',rank,'] ', kx
!!$ write(*,'(a,i5,a,16f10.4)') 'ky[',rank,'] ', ky
!!$
@@ -230,11 +230,11 @@ contains
!> forward transform - The result is saved in local buffers
- !! @param input data
+ !! @param input data
subroutine r2c_2d(inputData)
real(mk),dimension(:,:) :: inputData
-
+
integer(C_INTPTR_T) :: i, j
! init
do j = 1, local_resolution(c_Y)
@@ -253,7 +253,7 @@ contains
!!$ do i = 1, fft_resolution(c_Y)
!!$ write(*,'(a,i5,a,16f10.4)') 'aaaa[',rank,'] ', dataout1(i,1:local_resolution(c_X))
!!$ end do
-
+
end subroutine r2c_2d
!> Backward transform
@@ -280,7 +280,7 @@ contains
!!$ end do
!!$ write(*,'(a,i5,a)') '[',rank,'] ==============================='
-
+
end subroutine c2r_2d
!> Backward transform for scalar field
@@ -305,7 +305,7 @@ contains
!!$ end do
!!$ write(*,'(a,i5,a)') '[',rank,'] ==============================='
-
+
end subroutine c2r_scalar_2d
@@ -316,25 +316,25 @@ contains
!> Computation of frequencies coeff, over the distributed direction in the real/complex case
!> @param lengths size of the domain
subroutine computeKx(length)
-
+
real(mk),intent(in) :: length
-
+
!! Local loops indices
integer(C_INTPTR_T) :: i
!! Compute filter coefficients
allocate(kx(local_resolution(c_X)))
-
+
do i = local_offset(c_X)+1,local_offset(c_X)+local_resolution(c_X) !! global index
kx(i-local_offset(c_X)) = 2.*pi/length*(i-1)
end do
-
+
end subroutine computeKx
-
+
!> Computation of frequencies coeff, over the distributed direction in the complex/complex case
!> @param lengths size of the domain
subroutine computeKxC(length)
-
+
real(mk),intent(in) :: length
!! Local loops indices
@@ -342,17 +342,17 @@ contains
!! Compute filter coefficients
allocate(kx(local_resolution(c_X)))
-
+
!! x frequencies (distributed over proc)
!! If we deal with positive frequencies only
if(local_offset(c_X)+local_resolution(c_X) <= fft_resolution(c_X)/2+1 ) then
do i = 1,local_resolution(c_X)
kx(i) = 2.*pi/length*(local_offset(c_X)+i-1)
end do
-
+
else
!! else if we deal with negative frequencies only
- if(local_offset(c_X)+1 > fft_resolution(c_X)/2+1 ) then
+ if(local_offset(c_X)+1 > fft_resolution(c_X)/2+1 ) then
do i = 1,local_resolution(c_X)
kx(i) = 2.*pi/length*(local_offset(c_X)+i-1-fft_resolution(c_X))
end do
@@ -366,29 +366,29 @@ contains
end do
end if
end if
-
+
end subroutine computeKxC
!> Computation of frequencies coeff, over non-distributed direction(s)
!> @param lengths size of the domain
subroutine computeKy(length)
real(mk), intent(in) :: length
-
+
!! Local loops indices
integer(C_INTPTR_T) :: i
allocate(ky(fft_resolution(c_Y)))
-
+
do i = 1, fft_resolution(c_Y)/2+1
ky(i) = 2.*pi/length*(i-1)
end do
do i = fft_resolution(c_Y)/2+2,fft_resolution(c_Y)
ky(i) = 2.*pi/length*(i-fft_resolution(c_Y)-1)
end do
-
+
end subroutine computeKy
subroutine filter_poisson_2d()
-
+
integer(C_INTPTR_T) :: i, j
complex(C_DOUBLE_COMPLEX) :: coeff
if(local_offset(c_X)==0) then
@@ -426,28 +426,28 @@ contains
!!$ do i = 1,local_resolution(c_X)
!!$ write(*,'(a,i5,a,16f10.4)') 'xx[',rank,'] ', dataout1(1:fft_resolution(c_Y),i)
!!$ end do
-!!$
+!!$
!!$ write(*,'(a,i5,a)') '[',rank,'] ==============================='
!!$ do i = 1,local_resolution(c_X)
!!$ write(*,'(a,i5,a,16f10.4)') 'yy[',rank,'] ', dataout2(1:fft_resolution(c_Y),i)
!!$ end do
!!$ write(*,'(a,i5,a)') '[',rank,'] ==============================='
-
+
end subroutine filter_poisson_2d
-
+
subroutine filter_diffusion_2d(nudt)
-
+
real(C_DOUBLE), intent(in) :: nudt
integer(C_INTPTR_T) :: i, j
complex(C_DOUBLE_COMPLEX) :: coeff
-
+
do i = 1,local_resolution(c_X)
do j = 1, fft_resolution(c_Y)
coeff = 1./(1. + nudt * (kx(i)**2+ky(j)**2))
dataout1(j,i) = coeff*dataout1(j,i)
end do
end do
-
+
end subroutine filter_diffusion_2d
!> Clean fftw context (free memory, plans ...)
@@ -462,7 +462,7 @@ contains
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
@@ -472,11 +472,11 @@ contains
write(*,'(a,i5,a,2i12)') '[',rank,'] size of kx,y,z vectors (number of elements):', &
size(kx), size(ky)
write(*,'(a,i5,a,4i5)') '[',rank,'] local resolution and offset :', local_resolution, local_offset
- memoryAllocated = 2*memoryAllocated + real(sizeof(kx) + sizeof(ky))*1e-6
- write(*,'(a,i5,a,f10.2)') '[',rank,'] Total memory used for fftw buffers (MB):', memoryAllocated
+ memoryAllocated = 2*memoryAllocated + real(sizeof(kx) + sizeof(ky), mk)*1e-6
+ write(*,'(a,i5,a,f10.2)') '[',rank,'] Total memory used for fftw buffers (MB):', memoryAllocated
end subroutine fft2d_diagnostics
-
+
!> Get local size of input and output field on fftw topology
!! @param datashape local shape of the input field for the fftw process
!! @param offset index of the first component of the local field (in each dir) in the global set of indices
@@ -506,7 +506,7 @@ contains
howmany = 1
if(rank==0) open(unit=21,file=filename,form="formatted")
- print *, "ooOOOOOOOOOOOOO", FFTW_MPI_DEFAULT_BLOCK
+ print *, "ooOOOOOOOOOOOOO", FFTW_MPI_DEFAULT_BLOCK
allocate(fft_resolution(2))
fft_resolution(:) = resolution(:)-1
@@ -517,14 +517,14 @@ contains
alloc_local = fftw_mpi_local_size_many_transposed(2,n,howmany,FFTW_MPI_DEFAULT_BLOCK,&
FFTW_MPI_DEFAULT_BLOCK,MPI_COMM_WORLD,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)
! 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)])
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)
diff --git a/HySoP/src/fftw/fft3d.f90 b/HySoP/src/fftw/fft3d.f90
index 67093fc82..ea581db0e 100755
--- a/HySoP/src/fftw/fft3d.f90
+++ b/HySoP/src/fftw/fft3d.f90
@@ -7,48 +7,48 @@
!! input fields which are real. The names of these routines contain "r2c".
!! \li 3 - fftw routines for the "real to complex" case : solves the problem for
!! input fields which are real and using the "many" interface of the fftw.
-!! It means that transforms are applied to the 3 input fields at the same time.
-!! Names of these routines contain "many".
-!!
-!! Obviously, all the above cases should lead to the same results. By default
+!! It means that transforms are applied to the 3 input fields at the same time.
+!! Names of these routines contain "many".
+!!
+!! Obviously, all the above cases should lead to the same results. By default
!! case 2 must be chosen (if input is real). Case 1 and 3 are more or less
-!! dedicated to tests and validation.
+!! dedicated to tests and validation.
module fft3d
use client_data
- use mpi
+ use mpi
implicit none
include 'fftw3-mpi.f03'
- private
+ private
public :: init_r2c_3d,init_c2c_3d,r2c_3d,r2c_scalar_3d,c2c_3d,c2r_3d,c2r_scalar_3d,cleanFFTW_3d,&
getParamatersTopologyFFTW3d,filter_poisson_3d,filter_curl_diffusion_3d, &
init_r2c_3d_many, r2c_3d_many, c2r_3d_many, filter_diffusion_3d_many,&
filter_poisson_3d_many, filter_diffusion_3d, filter_curl_3d, filter_projection_om_3d,&
filter_multires_om_3d, filter_pressure_3d
-
+
!> plan for fftw "c2c" forward or r2c transform
type(C_PTR) :: plan_forward1, plan_forward2, plan_forward3
!> plan for fftw "c2c" backward or c2r transform
type(C_PTR) :: plan_backward1, plan_backward2, plan_backward3
!> memory buffer for fftw (input and output buffer will point to this location)
type(C_PTR) :: cbuffer1
- !> second memory buffer for fftw
+ !> second memory buffer for fftw
type(C_PTR) :: cbuffer2
- !> third memory buffer for fftw
+ !> third memory buffer for fftw
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_DOUBLE_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_DOUBLE), 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()
- !> Field (complex values) for fftw (forward) output
- complex(C_DOUBLE_COMPLEX), pointer :: dataout1(:,:,:)=>NULL() ,dataout2(:,:,:)=>NULL() ,dataout3(:,:,:)=>NULL()
+ real(C_DOUBLE), pointer :: rdatain_many(:,:,:,:)=>NULL()
+ !> Field (complex values) for fftw (forward) output
+ complex(C_DOUBLE_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_DOUBLE_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
@@ -72,7 +72,7 @@ module fft3d
contains
!========================================================================
- ! Complex to complex transforms
+ ! Complex to complex transforms
!========================================================================
!> Initialisation of the fftw context for complex to complex transforms (forward and backward)
@@ -82,29 +82,29 @@ contains
integer, dimension(3), intent(in) :: resolution
real(mk),dimension(3), intent(in) :: lengths
-
+
!! Size of the local memory required for fftw (cbuffer)
integer(C_INTPTR_T) :: alloc_local
-
+
if(is3DUpToDate) return
! init fftw mpi context
call fftw_mpi_init()
-
+
if(rank==0) open(unit=21,file=filename,form="formatted")
-
+
! set fft resolution
allocate(fft_resolution(3))
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),MPI_COMM_WORLD,&
local_resolution(c_Z),local_offset(c_Z),local_resolution(c_Y),local_offset(c_Y));
-
+
! Set a default value for c_X components.
local_offset(c_X) = 0
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)
! link datain and dataout to cbuffer, setting the right dimensions for each
@@ -136,7 +136,7 @@ contains
MPI_COMM_WORLD,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,&
MPI_COMM_WORLD,FFTW_BACKWARD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_IN))
-
+
call computeKx(lengths(c_X))
call computeKy(lengths(c_Y))
call computeKz(lengths(c_Z))
@@ -160,7 +160,7 @@ contains
subroutine c2c_3d(omega_x,omega_y,omega_z,velocity_x,velocity_y,velocity_z)
complex(mk),dimension(:,:,:) :: velocity_x,velocity_y,velocity_z
complex(mk),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
-
+
integer(C_INTPTR_T) :: i,j,k
! Copy input data into the fftw buffer
do k = 1, local_resolution(c_Z)
@@ -197,9 +197,9 @@ contains
end subroutine c2c_3d
!========================================================================
- ! Real to complex transforms
+ ! Real to complex transforms
!========================================================================
-
+
!> Initialisation of the fftw context for real to complex transforms (forward and backward)
!! @param[in] resolution global domain resolution
!! @param[in] lengths width of each side of the domain
@@ -214,16 +214,16 @@ contains
! init fftw mpi context
call fftw_mpi_init()
-
+
if(rank==0) open(unit=21,file=filename,form="formatted")
allocate(fft_resolution(3))
fft_resolution(:) = resolution(:)-1
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,&
MPI_COMM_WORLD,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.
local_offset(c_X) = 0
local_resolution(c_X) = halfLength
@@ -240,7 +240,7 @@ contains
call c_f_pointer(cbuffer2, dataout2, [halfLength, fft_resolution(c_Z), local_resolution(c_Y)])
call c_f_pointer(cbuffer3, rdatain3, [2*halfLength,fft_resolution(c_Y),local_resolution(c_Z)])
call c_f_pointer(cbuffer3, dataout3, [halfLength, fft_resolution(c_Z), local_resolution(c_Y)])
-
+
rdatain1 = 0.0
rdatain2 = 0.0
rdatain3 = 0.0
@@ -276,15 +276,15 @@ contains
!! @param[in] omega_z 3d scalar field, z-component of the input vector field
!! @param[in] ghosts, number of points in the ghost layer of input fields.
subroutine r2c_3d(omega_x,omega_y,omega_z, ghosts)
-
+
real(mk),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
integer, dimension(3), intent(in) :: ghosts
real(mk) :: start
integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
! ig, jg, kg are used to take into account
- ! ghost points in input fields
-
+ ! ghost points in input fields
+
! init
do k =1, local_resolution(c_Z)
kg = k + ghosts(c_Z)
@@ -307,7 +307,7 @@ contains
!!print *, "r2c time = ", MPI_WTIME() - start
end subroutine r2c_3d
-
+
!> Backward fftw transform
!! @param[in,out] velocity_x 3d scalar field, x-component of the output vector field
!! @param[in,out] velocity_y 3d scalar field, y-component of the output vector field
@@ -336,22 +336,22 @@ contains
end do
end do
end do
-
+
end subroutine c2r_3d
!> forward transform - The result is saved in a local buffer
!! @param[in] omega 3d scalar field, x-component of the input vector field
!! @param[in] ghosts, number of points in the ghost layer of input field.
subroutine r2c_scalar_3d(omega, ghosts)
-
+
real(mk),dimension(:,:,:),intent(in) :: omega
integer, dimension(3), intent(in) :: ghosts
real(mk) :: start
integer(C_INTPTR_T) :: i,j,k, ig, jg, kg
! ig, jg, kg are used to take into account
- ! ghost points in input fields
-
+ ! ghost points in input fields
+
! init
do k =1, local_resolution(c_Z)
kg = k + ghosts(c_Z)
@@ -370,7 +370,7 @@ contains
!!print *, "r2c time = ", MPI_WTIME() - start
end subroutine r2c_scalar_3d
-
+
!> Backward fftw transform
!! @param[in,out] omega 3d scalar field
!! @param[in] ghosts, number of points in the ghost layer of in/out omega scalar field.
@@ -393,13 +393,13 @@ contains
end do
end do
end do
-
+
end subroutine c2r_scalar_3d
!========================================================================
! Real to complex transforms based on "many" fftw routines
!========================================================================
-
+
!> Initialisation of the fftw context for real to complex transforms (forward and backward)
!! @param[in] resolution global domain resolution
!! @param[in] lengths width of each side of the domain
@@ -425,21 +425,21 @@ contains
! 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,&
MPI_COMM_WORLD,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.
local_offset(c_X) = 0
local_resolution(c_X) = halfLength
! allocate local buffer (used to save datain/dataout ==> in-place transform!!)
cbuffer1 = fftw_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)])
call c_f_pointer(cbuffer1, dataout_many, [howmany,halfLength, fft_resolution(c_Z), local_resolution(c_Y)])
-
+
! 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, &
MPI_COMM_WORLD,ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))
plan_backward1 = fftw_mpi_plan_many_dft_c2r(3,n,howmany,blocksize,blocksize, dataout_many, rdatain_many, &
@@ -460,13 +460,13 @@ contains
!! @param[in] omega_x 3d scalar field, x-component of the input vector field
!! @param[in] omega_y 3d scalar field, y-component of the input vector field
!! @param[in] omega_z 3d scalar field, z-component of the input vector field
- !! @param input data
+ !! @param input data
subroutine r2c_3d_many(omega_x,omega_y,omega_z)
real(mk),dimension(:,:,:),intent(in) :: omega_x,omega_y,omega_z
real(mk) :: start
integer(C_INTPTR_T) :: i,j,k
-
+
! init
do k =1, local_resolution(c_Z)
do j = 1, fft_resolution(c_Y)
@@ -477,14 +477,14 @@ contains
end do
end do
end do
-
+
! compute transform (as many times as desired)
start = MPI_WTIME()
call fftw_mpi_execute_dft_r2c(plan_forward1, rdatain_many, dataout_many)
!! print *, "r2c time = ", MPI_WTIME() - start
-
+
end subroutine r2c_3d_many
-
+
!> Backward fftw transform
!! @param[in,out] velocity_x 3d scalar field, x-component of the output vector field
!! @param[in,out] velocity_y 3d scalar field, y-component of the output vector field
@@ -508,7 +508,7 @@ contains
end do
end subroutine c2r_3d_many
-
+
!========================================================================
! Common (r2c, c2c) subroutines
!========================================================================
@@ -516,9 +516,9 @@ contains
!> Computation of frequencies coeff, over the distributed direction in the real/complex case
!> @param lengths size of the domain
subroutine computeKx(length)
-
+
real(mk),intent(in) :: length
-
+
!! Local loops indices
integer(C_INTPTR_T) :: i
@@ -532,12 +532,12 @@ contains
kx(i) = 2.*pi/length*(i-fft_resolution(c_X)-1)
end do
end subroutine computeKx
-
+
!> Computation of frequencies coeff, over distributed direction(s)
!> @param lengths size of the domain
subroutine computeKy(length)
real(C_DOUBLE), intent(in) :: length
-
+
!! Local loops indices
integer(C_INTPTR_T) :: i
allocate(ky(local_resolution(c_Y)))
@@ -550,7 +550,7 @@ contains
end do
else
!! else if we deal with negative frequencies only
- if(local_offset(c_Y)+1 > fft_resolution(c_Y)/2+1 ) then
+ if(local_offset(c_Y)+1 > fft_resolution(c_Y)/2+1 ) then
do i = 1,local_resolution(c_Y)
ky(i) = 2.*pi/length*(local_offset(c_Y)+i-1-fft_resolution(c_Y))
end do
@@ -564,14 +564,14 @@ contains
end do
end if
end if
-
+
end subroutine computeKy
!> Computation of frequencies coeff, over non-distributed direction(s)
!> @param length size of the domain
subroutine computeKz(length)
real(mk),intent(in) :: length
-
+
!! Local loops indices
integer(C_INTPTR_T) :: i
allocate(kz(fft_resolution(c_Z)))
@@ -583,15 +583,15 @@ contains
end do
end subroutine computeKz
-
+
!> Solve Poisson problem in the Fourier space :
!! \f{eqnarray*} \Delta \psi &=& - \omega \\ v = \nabla\times\psi \f}
subroutine filter_poisson_3d()
-
+
integer(C_INTPTR_T) :: i,j,k
complex(C_DOUBLE_COMPLEX) :: coeff
complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2
-
+
! Set first coeff (check for "all freq = 0" case)
if(local_offset(c_Y) == 0) then
dataout1(1,1,1) = 0.0
@@ -599,12 +599,12 @@ contains
dataout3(1,1,1) = 0.0
else
coeff = Icmplx/(ky(1)**2)
- buffer1 = dataout1(1,1,1)
+ buffer1 = dataout1(1,1,1)
dataout1(1,1,1) = coeff*ky(1)*dataout3(1,1,1)
dataout2(1,1,1) = 0.0
dataout3(1,1,1) = -coeff*ky(1)*buffer1
endif
-
+
! !! mind the transpose -> index inversion between y and z
do i = 2, local_resolution(c_X)
coeff = Icmplx/(kx(i)**2+ky(1)**2)
@@ -619,14 +619,14 @@ contains
do k = 2, fft_resolution(c_Z)
do i = 1, local_resolution(c_X)
coeff = Icmplx/(kx(i)**2+ky(1)**2+kz(k)**2)
- buffer1 = dataout1(i,k,1)
+ buffer1 = dataout1(i,k,1)
buffer2 = dataout2(i,k,1)
- dataout1(i,k,1) = coeff*(ky(1)*dataout3(i,k,1)-kz(k)*dataout2(i,k,1))
+ dataout1(i,k,1) = coeff*(ky(1)*dataout3(i,k,1)-kz(k)*dataout2(i,k,1))
dataout2(i,k,1) = coeff*(kz(k)*buffer1-kx(i)*dataout3(i,k,1))
dataout3(i,k,1) = coeff*(kx(i)*buffer2-ky(1)*buffer1)
end do
end do
-
+
! !! mind the transpose -> index inversion between y and z
do j = 2,local_resolution(c_Y)
do k = 1, fft_resolution(c_Z)
@@ -634,15 +634,15 @@ contains
coeff = Icmplx/(kx(i)**2+ky(j)**2+kz(k)**2)
buffer1 = dataout1(i,k,j)
buffer2 = dataout2(i,k,j)
- dataout1(i,k,j) = coeff*(ky(j)*dataout3(i,k,j)-kz(k)*dataout2(i,k,j))
+ dataout1(i,k,j) = coeff*(ky(j)*dataout3(i,k,j)-kz(k)*dataout2(i,k,j))
dataout2(i,k,j) = coeff*(kz(k)*buffer1-kx(i)*dataout3(i,k,j))
dataout3(i,k,j) = coeff*(kx(i)*buffer2-ky(j)*buffer1)
end do
end do
end do
-
+
end subroutine filter_poisson_3d
-
+
!> Solve diffusion problem in the Fourier space :
!! \f{eqnarray*} \omega &=& \nabla \times v \\ \frac{\partial \omega}{\partial t} &=& \nu \Delta \omega \f}
!! @param[in] nudt \f$ \nu\times dt\f$, diffusion coefficient times current time step
@@ -652,7 +652,7 @@ contains
integer(C_INTPTR_T) :: i,j,k
complex(C_DOUBLE_COMPLEX) :: coeff
complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2
-
+
!! mind the transpose -> index inversion between y and z
do j = 1,local_resolution(c_Y)
do k = 1, fft_resolution(c_Z)
@@ -660,7 +660,7 @@ contains
coeff = Icmplx/(1. + nudt * (kx(i)**2+ky(j)**2+kz(k)**2))
buffer1 = dataout1(i,k,j)
buffer2 = dataout2(i,k,j)
- dataout1(i,k,j) = coeff*(ky(j)*dataout3(i,k,j)-kz(k)*dataout2(i,k,j))
+ dataout1(i,k,j) = coeff*(ky(j)*dataout3(i,k,j)-kz(k)*dataout2(i,k,j))
dataout2(i,k,j) = coeff*(kz(k)*buffer1-kx(i)*dataout3(i,k,j))
dataout3(i,k,j) = coeff*(kx(i)*buffer2-ky(j)*buffer1)
end do
@@ -677,7 +677,7 @@ contains
real(C_DOUBLE), intent(in) :: nudt
integer(C_INTPTR_T) :: i,j,k
complex(C_DOUBLE_COMPLEX) :: coeff
-
+
!! mind the transpose -> index inversion between y and z
do j = 1,local_resolution(c_Y)
do k = 1, fft_resolution(c_Z)
@@ -699,7 +699,7 @@ contains
integer(C_INTPTR_T) :: i,j,k
complex(C_DOUBLE_COMPLEX) :: coeff
complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2
-
+
!! mind the transpose -> index inversion between y and z
do j = 1,local_resolution(c_Y)
do k = 1, fft_resolution(c_Z)
@@ -719,11 +719,11 @@ contains
!> Perform solenoidal projection to ensure divergence free vorticity field
!! \f{eqnarray*} \omega ' &=& \omega - \nabla\pi \f}
subroutine filter_projection_om_3d()
-
+
integer(C_INTPTR_T) :: i,j,k
complex(C_DOUBLE_COMPLEX) :: coeff
complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2,buffer3
-
+
! Set first coeff (check for "all freq = 0" case)
if(local_offset(c_Y) == 0) then
dataout1(1,1,1) = dataout1(1,1,1)
@@ -736,7 +736,7 @@ contains
dataout2(1,1,1) = dataout2(1,1,1) - coeff*ky(1)*(ky(1)*buffer2)
dataout3(1,1,1) = dataout3(1,1,1)
endif
-
+
! !! mind the transpose -> index inversion between y and z
do i = 2, local_resolution(c_X)
coeff = 1./(kx(i)**2+ky(1)**2+kz(1)**2)
@@ -760,7 +760,7 @@ contains
dataout3(i,k,1) = dataout3(i,k,1) - coeff*kz(k)*(kx(i)*buffer1+ky(1)*buffer2+kz(k)*buffer3)
end do
end do
-
+
! !! mind the transpose -> index inversion between y and z
do j = 2,local_resolution(c_Y)
do k = 1, fft_resolution(c_Z)
@@ -775,7 +775,7 @@ contains
end do
end do
end do
-
+
end subroutine filter_projection_om_3d
!> Projects vorticity values from fine to coarse grid :
@@ -806,12 +806,12 @@ contains
end subroutine filter_multires_om_3d
- !> Solve the Poisson problem allowing to recover
+ !> Solve the Poisson problem allowing to recover
!! pressure from velocity in the Fourier space
subroutine filter_pressure_3d()
integer(C_INTPTR_T) :: i,j,k
complex(C_DOUBLE_COMPLEX) :: coeff
-
+
! Set first coeff (check for "all freq = 0" case)
if(local_offset(c_Y) == 0) then
dataout1(1,1,1) = 0.0
@@ -819,7 +819,7 @@ contains
coeff = -1./(ky(1)**2)
dataout1(1,1,1) = coeff*dataout1(1,1,1)
endif
-
+
! !! mind the transpose -> index inversion between y and z
do i = 2, local_resolution(c_X)
coeff = -1./(kx(i)**2+ky(1)**2)
@@ -833,7 +833,7 @@ contains
dataout1(i,k,1) = coeff*dataout1(i,k,1)
end do
end do
-
+
! !! mind the transpose -> index inversion between y and z
do j = 2,local_resolution(c_Y)
do k = 1, fft_resolution(c_Z)
@@ -845,26 +845,26 @@ contains
end do
end subroutine filter_pressure_3d
-
+
!> Solve Poisson problem in the Fourier space :
!! \f{eqnarray*} \Delta \psi &=& - \omega \\ v &=& \nabla\times\psi \f}
subroutine filter_poisson_3d_many()
-
+
integer(C_INTPTR_T) :: i,j,k
complex(C_DOUBLE_COMPLEX) :: coeff
complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2
-
+
! Set first coeff (check for "all freq = 0" case)
if(local_offset(c_Y) == 0) then
dataout_many(:,1,1,1) = 0.0
else
coeff = Icmplx/(ky(1)**2)
- buffer1 = dataout_many(1,1,1,1)
+ buffer1 = dataout_many(1,1,1,1)
dataout_many(1,1,1,1) = coeff*ky(1)*dataout_many(3,1,1,1)
dataout_many(2,1,1,1) = 0.0
dataout_many(3,1,1,1) = -coeff*ky(1)*buffer1
endif
-
+
! !! mind the transpose -> index inversion between y and z
do i = 2, local_resolution(c_X)
coeff = Icmplx/(kx(i)**2+ky(1)**2)
@@ -879,14 +879,14 @@ contains
do k = 2, fft_resolution(c_Z)
do i = 1, local_resolution(c_X)
coeff = Icmplx/(kx(i)**2+ky(1)**2+kz(k)**2)
- buffer1 = dataout_many(1,i,k,1)
+ buffer1 = dataout_many(1,i,k,1)
buffer2 = dataout_many(2,i,k,1)
- dataout_many(1,i,k,1) = coeff*(ky(1)*dataout_many(3,i,k,1)-kz(k)*dataout_many(2,i,k,1))
+ dataout_many(1,i,k,1) = coeff*(ky(1)*dataout_many(3,i,k,1)-kz(k)*dataout_many(2,i,k,1))
dataout_many(2,i,k,1) = coeff*(kz(k)*buffer1-kx(i)*dataout_many(3,i,k,1))
dataout_many(3,i,k,1) = coeff*(kx(i)*buffer2-ky(1)*buffer1)
end do
end do
-
+
! !! mind the transpose -> index inversion between y and z
do j = 2,local_resolution(c_Y)
do k = 1, fft_resolution(c_Z)
@@ -894,15 +894,15 @@ contains
coeff = Icmplx/(kx(i)**2+ky(j)**2+kz(k)**2)
buffer1 = dataout_many(1,i,k,j)
buffer2 = dataout_many(2,i,k,j)
- dataout_many(1,i,k,j) = coeff*(ky(j)*dataout_many(3,i,k,j)-kz(k)*dataout_many(2,i,k,j))
+ dataout_many(1,i,k,j) = coeff*(ky(j)*dataout_many(3,i,k,j)-kz(k)*dataout_many(2,i,k,j))
dataout_many(2,i,k,j) = coeff*(kz(k)*buffer1-kx(i)*dataout_many(3,i,k,j))
dataout_many(3,i,k,j) = coeff*(kx(i)*buffer2-ky(j)*buffer1)
end do
end do
end do
-
+
end subroutine filter_poisson_3d_many
-
+
!> Solve diffusion problem in the Fourier space :
!! \f{eqnarray*} \omega &=& \nabla \times v \\ \frac{\partial \omega}{\partial t} &=& \nu \Delta \omega \f}
!! @param[in] nudt \f$ \nu\times dt\f$, diffusion coefficient times current time step
@@ -912,7 +912,7 @@ contains
integer(C_INTPTR_T) :: i,j,k
complex(C_DOUBLE_COMPLEX) :: coeff
complex(C_DOUBLE_COMPLEX) :: buffer1,buffer2
-
+
!! mind the transpose -> index inversion between y and z
do j = 1,local_resolution(c_Y)
do k = 1, fft_resolution(c_Z)
@@ -920,7 +920,7 @@ contains
coeff = Icmplx/(1. + nudt * kx(i)**2+ky(j)**2+kz(k)**2)
buffer1 = dataout_many(1,i,k,j)
buffer2 = dataout_many(2,i,k,j)
- dataout_many(1,i,k,j) = coeff*(ky(j)*dataout_many(3,i,k,j)-kz(k)*dataout_many(2,i,k,j))
+ dataout_many(1,i,k,j) = coeff*(ky(j)*dataout_many(3,i,k,j)-kz(k)*dataout_many(2,i,k,j))
dataout_many(2,i,k,j) = coeff*(kz(k)*buffer1-kx(i)*dataout_many(3,i,k,j))
dataout_many(3,i,k,j) = coeff*(kx(i)*buffer2-ky(j)*buffer1)
end do
@@ -951,7 +951,7 @@ contains
!> some information about memory alloc, arrays sizes and so on
subroutine fft3d_diagnostics(nbelem,howmany)
integer(C_INTPTR_T), intent(in) :: nbelem
- ! number of buffers used for fftw
+ ! number of buffers used for fftw
integer, optional,intent(in) :: howmany
complex(C_DOUBLE_COMPLEX) :: memoryAllocated
@@ -967,8 +967,8 @@ contains
write(*,'(a,i5,a,3i12)') '[',rank,'] size of kx,y,z vectors (number of elements):', &
size(kx), size(ky),size(kz)
write(*,'(a,i5,a,6i5)') '[',rank,'] local resolution and offset :', local_resolution, local_offset
- memoryAllocated = nbFields*memoryAllocated + real(sizeof(kx) + sizeof(ky) + sizeof(kz))*1e-6
- write(*,'(a,i5,a,f10.2)') '[',rank,'] Total memory used for fftw buffers (MB):', memoryAllocated
+ memoryAllocated = nbFields*memoryAllocated + real(sizeof(kx) + sizeof(ky) + sizeof(kz), mk)*1e-6
+ write(*,'(a,i5,a,f10.2)') '[',rank,'] Total memory used for fftw buffers (MB):', memoryAllocated
end subroutine fft3d_diagnostics
@@ -981,8 +981,8 @@ contains
integer(C_INTPTR_T) :: zero = 0
datashape = (/fft_resolution(c_X), fft_resolution(c_Y), local_resolution(c_Z)/)
offset = (/zero, zero, local_offset(c_Z)/)
-
+
end subroutine getParamatersTopologyFFTW3d
-
+
end module fft3d
diff --git a/HySoP/src/scalesInterface/particles/advec.f90 b/HySoP/src/scalesInterface/particles/advec.f90
index 10303ff06..a34d30614 100644
--- a/HySoP/src/scalesInterface/particles/advec.f90
+++ b/HySoP/src/scalesInterface/particles/advec.f90
@@ -119,12 +119,12 @@ subroutine advec_init(order, stab_coeff, verbosity, dim_split)
case('classic')
advec_step => advec_step_Torder1
! Compute stability coefficient
- if (present(stab_coeff)) stab_coeff = 1.0/(2.0*dble(bl_size))
+ if (present(stab_coeff)) stab_coeff = 1.0/(2.0*real(bl_size, WP))
case default ! Strang
advec_step => advec_step_Torder2
! Compute stability coefficient - as each dimension is solved in
! dt/2, stab_coef is 2 times bigger
- if (present(stab_coeff)) stab_coeff = 1.0/(dble(bl_size))
+ if (present(stab_coeff)) stab_coeff = 1.0/(real(bl_size, WP))
end select
! Call the right remeshing formula
diff --git a/HySoP/src/scalesInterface/particles/advecX.f90 b/HySoP/src/scalesInterface/particles/advecX.f90
index 33714e663..9cb350e43 100644
--- a/HySoP/src/scalesInterface/particles/advecX.f90
+++ b/HySoP/src/scalesInterface/particles/advecX.f90
@@ -167,7 +167,7 @@ subroutine advecX_remesh_in_buffer_lambda(gs, j, k, ind_min, p_pos_adim, bl_type
! -- Allocate remeshX_pter --
allocate(remeshX_pter(send_j_min:send_j_max))
do ind = send_j_min, send_j_max
- proc_gap = floor(real(ind-1)/mesh_sc%N_proc(direction)) - (ind_min-1)
+ proc_gap = floor(real(ind-1, WP)/mesh_sc%N_proc(direction)) - (ind_min-1)
remeshX_pter(ind)%pter => buffer(pos_in_buffer(proc_gap))
pos_in_buffer(proc_gap) = pos_in_buffer(proc_gap) + 1
end do
@@ -242,7 +242,7 @@ subroutine advecX_remesh_in_buffer_limit_lambda(gs, j, k, ind_min, p_pos_adim, b
! -- Allocate remeshX_pter --
allocate(remeshX_pter(send_j_min:send_j_max))
do ind = send_j_min, send_j_max
- proc_gap = floor(real(ind-1)/mesh_sc%N_proc(direction)) - (ind_min-1)
+ proc_gap = floor(real(ind-1, WP)/mesh_sc%N_proc(direction)) - (ind_min-1)
remeshX_pter(ind)%pter => buffer(pos_in_buffer(proc_gap))
pos_in_buffer(proc_gap) = pos_in_buffer(proc_gap) + 1
end do
@@ -298,8 +298,8 @@ subroutine advecX_remesh_in_buffer_Mprime(gs, j, k, ind_min, p_pos_adim, send_mi
! sorted by receivers
integer :: i1, i2 ! indice of a line into the group
integer :: ind ! indice of the current particle inside the current line.
- real(WP), dimension(mesh_sc%N_proc(direction)) :: pos_translat ! translation of p_pos_adim as array indice
- ! are now starting from 1 and not ind_min
+ !! real(WP), dimension(mesh_sc%N_proc(direction)) :: pos_translat ! translation of p_pos_adim as array indice
+ !! ! are now starting from 1 and not ind_min
! ===== Remeshing into the buffer by using pointer array =====
@@ -309,16 +309,18 @@ subroutine advecX_remesh_in_buffer_Mprime(gs, j, k, ind_min, p_pos_adim, send_mi
! -- Allocate remeshX_pter --
allocate(remeshX_pter(send_min(i1,i2):send_max(i1,i2)))
do ind = send_min(i1,i2), send_max(i1,i2)
- proc_gap = floor(real(ind-1)/mesh_sc%N_proc(direction)) - (ind_min-1)
+ proc_gap = floor(real(ind-1, WP)/mesh_sc%N_proc(direction)) - (ind_min-1)
remeshX_pter(ind)%pter => buffer(pos_in_buffer(proc_gap))
pos_in_buffer(proc_gap) = pos_in_buffer(proc_gap) + 1
end do
- pos_translat = p_pos_adim(:,i1,i2) - send_min(i1,i2) + 1
+ !! pos_translat = p_pos_adim(:,i1,i2) - send_min(i1,i2) + 1
+ !! Index translation is performed in the AC_remesh_Mprime_pter subroutine on the
+ !! integer adimensionned particle position instead of here on the float position
! -- Remesh the particles in the buffer --
do ind = 1, mesh_sc%N_proc(direction)
- call AC_remesh_Mprime_pter(pos_translat(ind), scalar(ind,j+i1-1,k+i2-1), remeshX_pter)
+ call AC_remesh_Mprime_pter(p_pos_adim(ind,i1,i2), 1-send_min(i1,i2), scalar(ind,j+i1-1,k+i2-1), remeshX_pter)
end do
deallocate(remeshX_pter)
diff --git a/HySoP/src/scalesInterface/particles/advecY.f90 b/HySoP/src/scalesInterface/particles/advecY.f90
index 1bd653921..c7c72bbcd 100644
--- a/HySoP/src/scalesInterface/particles/advecY.f90
+++ b/HySoP/src/scalesInterface/particles/advecY.f90
@@ -121,7 +121,7 @@ subroutine advecY_remesh_in_buffer_lambda(gs, i, k, ind_min, p_pos_adim, bl_type
! -- Allocate remeshY_pter --
allocate(remeshY_pter(send_j_min:send_j_max))
do ind = send_j_min, send_j_max
- proc_gap = floor(real(ind-1)/mesh_sc%N_proc(direction)) - (ind_min-1)
+ proc_gap = floor(real(ind-1, WP)/mesh_sc%N_proc(direction)) - (ind_min-1)
remeshY_pter(ind)%pter => buffer(pos_in_buffer(proc_gap))
pos_in_buffer(proc_gap) = pos_in_buffer(proc_gap) + 1
end do
@@ -197,7 +197,7 @@ subroutine advecY_remesh_in_buffer_limit_lambda(gs, i, k, ind_min, p_pos_adim, b
! -- Allocate remeshY_pter --
allocate(remeshY_pter(send_j_min:send_j_max))
do ind = send_j_min, send_j_max
- proc_gap = floor(real(ind-1)/mesh_sc%N_proc(direction)) - (ind_min-1)
+ proc_gap = floor(real(ind-1, WP)/mesh_sc%N_proc(direction)) - (ind_min-1)
remeshY_pter(ind)%pter => buffer(pos_in_buffer(proc_gap))
pos_in_buffer(proc_gap) = pos_in_buffer(proc_gap) + 1
end do
@@ -253,8 +253,8 @@ subroutine advecY_remesh_in_buffer_Mprime(gs, i, k, ind_min, p_pos_adim, send_mi
! sorted by receivers
integer :: i1, i2 ! indice of a line into the group
integer :: ind ! indice of the current particle inside the current line.
- real(WP), dimension(mesh_sc%N_proc(direction)) :: pos_translat ! translation of p_pos_adim as array indice
- ! are now starting from 1 and not ind_min
+ !! real(WP), dimension(mesh_sc%N_proc(direction)) :: pos_translat ! translation of p_pos_adim as array indice
+ !! ! are now starting from 1 and not ind_min
! ===== Remeshing into the buffer by using pointer array =====
@@ -264,16 +264,18 @@ subroutine advecY_remesh_in_buffer_Mprime(gs, i, k, ind_min, p_pos_adim, send_mi
! -- Allocate remeshX_pter --
allocate(remeshY_pter(send_min(i1,i2):send_max(i1,i2)))
do ind = send_min(i1,i2), send_max(i1,i2)
- proc_gap = floor(real(ind-1)/mesh_sc%N_proc(direction)) - (ind_min-1)
+ proc_gap = floor(real(ind-1, WP)/mesh_sc%N_proc(direction)) - (ind_min-1)
remeshY_pter(ind)%pter => buffer(pos_in_buffer(proc_gap))
pos_in_buffer(proc_gap) = pos_in_buffer(proc_gap) + 1
end do
- pos_translat = p_pos_adim(:,i1,i2) - send_min(i1,i2) + 1
+ !! pos_translat = p_pos_adim(:,i1,i2) - send_min(i1,i2) + 1
+ !! Index translation is performed in the AC_remesh_Mprime_pter subroutine on the
+ !! integer adimensionned particle position instead of here on the float position
! -- Remesh the particles in the buffer --
do ind = 1, mesh_sc%N_proc(direction)
- call AC_remesh_Mprime_pter(pos_translat(ind), scalar(i+i1-1,ind,k+i2-1), remeshY_pter)
+ call AC_remesh_Mprime_pter(p_pos_adim(ind,i1,i2), 1-send_min(i1,i2), scalar(i+i1-1,ind,k+i2-1), remeshY_pter)
end do
deallocate(remeshY_pter)
diff --git a/HySoP/src/scalesInterface/particles/advecZ.f90 b/HySoP/src/scalesInterface/particles/advecZ.f90
index 898c4abea..e5c5d6f17 100644
--- a/HySoP/src/scalesInterface/particles/advecZ.f90
+++ b/HySoP/src/scalesInterface/particles/advecZ.f90
@@ -119,7 +119,7 @@ subroutine advecZ_remesh_in_buffer_lambda(gs, i, j, ind_min, p_pos_adim, bl_type
! -- Allocate remeshX_pter --
allocate(remeshZ_pter(send_j_min:send_j_max))
do ind = send_j_min, send_j_max
- proc_gap = floor(real(ind-1)/mesh_sc%N_proc(direction)) - (ind_min-1)
+ proc_gap = floor(real(ind-1, WP)/mesh_sc%N_proc(direction)) - (ind_min-1)
remeshZ_pter(ind)%pter => buffer(pos_in_buffer(proc_gap))
pos_in_buffer(proc_gap) = pos_in_buffer(proc_gap) + 1
end do
@@ -194,7 +194,7 @@ subroutine advecZ_remesh_in_buffer_limit_lambda(gs, i, j, ind_min, p_pos_adim, b
! -- Allocate remeshX_pter --
allocate(remeshZ_pter(send_j_min:send_j_max))
do ind = send_j_min, send_j_max
- proc_gap = floor(real(ind-1)/mesh_sc%N_proc(direction)) - (ind_min-1)
+ proc_gap = floor(real(ind-1, WP)/mesh_sc%N_proc(direction)) - (ind_min-1)
remeshZ_pter(ind)%pter => buffer(pos_in_buffer(proc_gap))
pos_in_buffer(proc_gap) = pos_in_buffer(proc_gap) + 1
end do
@@ -252,8 +252,8 @@ subroutine advecZ_remesh_in_buffer_Mprime(gs, i, j, ind_min, p_pos_adim, send_mi
! sorted by receivers
integer :: i1, i2 ! indice of a line into the group
integer :: ind ! indice of the current particle inside the current line.
- real(WP), dimension(mesh_sc%N_proc(direction)) :: pos_translat ! translation of p_pos_adim as array indice
- ! are now starting from 1 and not ind_min
+ !! real(WP), dimension(mesh_sc%N_proc(direction)) :: pos_translat ! translation of p_pos_adim as array indice
+ !! ! are now starting from 1 and not ind_min
! ===== Remeshing into the buffer by using pointer array =====
@@ -263,16 +263,18 @@ subroutine advecZ_remesh_in_buffer_Mprime(gs, i, j, ind_min, p_pos_adim, send_mi
! -- Allocate remeshZ_pter --
allocate(remeshZ_pter(send_min(i1,i2):send_max(i1,i2)))
do ind = send_min(i1,i2), send_max(i1,i2)
- proc_gap = floor(real(ind-1)/mesh_sc%N_proc(direction)) - (ind_min-1)
+ proc_gap = floor(real(ind-1, WP)/mesh_sc%N_proc(direction)) - (ind_min-1)
remeshZ_pter(ind)%pter => buffer(pos_in_buffer(proc_gap))
pos_in_buffer(proc_gap) = pos_in_buffer(proc_gap) + 1
end do
- pos_translat = p_pos_adim(:,i1,i2) - send_min(i1,i2) + 1
+ !! pos_translat = p_pos_adim(:,i1,i2) - send_min(i1,i2) + 1
+ !! Index translation is performed in the AC_remesh_Mprime_pter subroutine on the
+ !! integer adimensionned particle position instead of here on the float position
! -- Remesh the particles in the buffer --
do ind = 1, mesh_sc%N_proc(direction)
- call AC_remesh_Mprime_pter(pos_translat(ind), scalar(i+i1-1,j+i2-1,ind), remeshZ_pter)
+ call AC_remesh_Mprime_pter(p_pos_adim(ind,i1,i2), 1-send_min(i1,i2), scalar(i+i1-1,j+i2-1,ind), remeshZ_pter)
end do
deallocate(remeshZ_pter)
diff --git a/HySoP/src/scalesInterface/particles/advec_common_interpol.F90 b/HySoP/src/scalesInterface/particles/advec_common_interpol.F90
index 0398c6270..defedc610 100644
--- a/HySoP/src/scalesInterface/particles/advec_common_interpol.F90
+++ b/HySoP/src/scalesInterface/particles/advec_common_interpol.F90
@@ -155,8 +155,8 @@ subroutine AC_interpol_lin(direction, gs, ind_group, V_comp, p_inter)
! ===== Exchange velocity field if needed =====
! It uses non blocking message to do the computations during the communication process
! -- What have I to communicate ? --
- rece_gap(:,:,1) = floor(real(rece_ind_min-1)/mesh_sc%N_proc(direction))
- rece_gap(:,:,2) = floor(real(rece_ind_max-1)/mesh_sc%N_proc(direction))
+ rece_gap(:,:,1) = floor(real(rece_ind_min-1, WP)/mesh_sc%N_proc(direction))
+ rece_gap(:,:,2) = floor(real(rece_ind_max-1, WP)/mesh_sc%N_proc(direction))
rece_gap_abs(1) = minval(rece_gap(:,:,1))
rece_gap_abs(2) = maxval(rece_gap(:,:,2))
max_size = 2 + gs(2)*(2+3*gs(1))
@@ -334,7 +334,7 @@ subroutine AC_interpol_lin(direction, gs, ind_group, V_comp, p_inter)
weight = p_inter(ind,i1,i2)-pos
#endif
! Vm = V(pos)
- proc_gap = floor(real(pos-1)/mesh_sc%N_proc(direction))
+ proc_gap = floor(real(pos-1, WP)/mesh_sc%N_proc(direction))
if (neighbors(direction,proc_gap) == D_rank(direction)) then
#ifndef BLOCKING_SEND_PLUS
Vm(ind,i1,i2)%pter => V_comp(pos-proc_gap*mesh_sc%N_proc(direction), i1,i2)
@@ -352,7 +352,7 @@ subroutine AC_interpol_lin(direction, gs, ind_group, V_comp, p_inter)
myself(1) = .false.
end if
! Vp = V(pos+1)
- gap = floor(real(pos+1-1)/mesh_sc%N_proc(direction))
+ gap = floor(real(pos+1-1, WP)/mesh_sc%N_proc(direction))
if (neighbors(direction,gap) == D_rank(direction)) then
#ifndef BLOCKING_SEND_PLUS
Vp(ind,i1,i2)%pter => V_comp(pos+1-gap*mesh_sc%N_proc(direction), i1,i2)
@@ -434,7 +434,7 @@ subroutine AC_interpol_lin(direction, gs, ind_group, V_comp, p_inter)
! -- When we reach the end of the sub-domain OR the end of the particle line --
if (pos>proc_end) then ! Changement of subdomain
! We have reach the next subdomain => update values
- proc_gap = floor(real(pos-1)/mesh_sc%N_proc(direction)) ! "proc_gap = proc_gap + 1" does not work if N_proc = 1 and pos-pos_old = 2.
+ proc_gap = floor(real(pos-1, WP)/mesh_sc%N_proc(direction)) ! "proc_gap = proc_gap + 1" does not work if N_proc = 1 and pos-pos_old = 2.
myself(1) = (neighbors(direction,proc_gap) == D_rank(direction)) ! For the same reason that line jsute above, we do not use "myself(1) = myself(2)"
proc_end = (proc_gap+1)*mesh_sc%N_proc(direction)
myself(2) = (neighbors(direction,proc_gap+1) == D_rank(direction))
@@ -539,7 +539,7 @@ subroutine AC_interpol_lin(direction, gs, ind_group, V_comp, p_inter)
! -- When we reach the end of the sub-domain OR the end of the particle line --
if (pos>proc_end) then ! Changement of subdomain
! We have reach the next subdomain => update values
- proc_gap = floor(real(pos-1)/mesh_sc%N_proc(direction)) ! "proc_gap = proc_gap + 1" does not work if N_proc = 1 and pos-pos_old = 2.
+ proc_gap = floor(real(pos-1, WP)/mesh_sc%N_proc(direction)) ! "proc_gap = proc_gap + 1" does not work if N_proc = 1 and pos-pos_old = 2.
myself(1) = (neighbors(direction,proc_gap) == D_rank(direction)) ! For the same reason that line jsute above, we do not use "myself(1) = myself(2)"
proc_end = (proc_gap+1)*mesh_sc%N_proc(direction)
myself(2) = (neighbors(direction,proc_gap+1) == D_rank(direction))
@@ -963,8 +963,8 @@ subroutine AC_interpol_plus(direction, gs, ind_group, id1, id2, V_coarse, p_V)
! ===== Exchange velocity field if needed =====
! It uses non blocking message to do the computations during the communication process
! -- What have I to communicate ? --
- rece_gap(:,:,1) = floor(real(rece_ind_min-1)/mesh_V%N_proc(direction))
- rece_gap(:,:,2) = floor(real(rece_ind_max-1)/mesh_V%N_proc(direction))
+ rece_gap(:,:,1) = floor(real(rece_ind_min-1, WP)/mesh_V%N_proc(direction))
+ rece_gap(:,:,2) = floor(real(rece_ind_max-1, WP)/mesh_V%N_proc(direction))
rece_gap_abs(1) = minval(rece_gap(:,:,1))
rece_gap_abs(2) = maxval(rece_gap(:,:,2))
max_size = 2 + gs(2)*(2+3*gs(1))
diff --git a/HySoP/src/scalesInterface/particles/advec_common_remesh.F90 b/HySoP/src/scalesInterface/particles/advec_common_remesh.F90
index 371284faf..83316eec6 100644
--- a/HySoP/src/scalesInterface/particles/advec_common_remesh.F90
+++ b/HySoP/src/scalesInterface/particles/advec_common_remesh.F90
@@ -813,8 +813,8 @@ subroutine AC_remesh_range(bl_type, p_pos_adim, direction, send_min, send_max, s
end where
! -- What have I to communicate ? --
- send_gap(:,:,1) = floor(real(send_min-1)/mesh_sc%N_proc(direction))
- send_gap(:,:,2) = floor(real(send_max-1)/mesh_sc%N_proc(direction))
+ send_gap(:,:,1) = floor(real(send_min-1, WP)/mesh_sc%N_proc(direction))
+ send_gap(:,:,2) = floor(real(send_max-1, WP)/mesh_sc%N_proc(direction))
send_gap_abs(1) = minval(send_gap(:,:,1))
send_gap_abs(2) = maxval(send_gap(:,:,2))
@@ -847,8 +847,8 @@ subroutine AC_remesh_range_notype(p_pos_adim, direction, send_min, send_max, sen
send_max = floor(p_pos_adim(mesh_sc%N_proc(direction),:,:))+remesh_stencil(2)
! -- What have I to communicate ? --
- send_gap(:,:,1) = floor(real(send_min-1)/mesh_sc%N_proc(direction))
- send_gap(:,:,2) = floor(real(send_max-1)/mesh_sc%N_proc(direction))
+ send_gap(:,:,1) = floor(real(send_min-1, WP)/mesh_sc%N_proc(direction))
+ send_gap(:,:,2) = floor(real(send_max-1, WP)/mesh_sc%N_proc(direction))
send_gap_abs(1) = minval(send_gap(:,:,1))
send_gap_abs(2) = maxval(send_gap(:,:,2))
diff --git a/HySoP/src/scalesInterface/particles/advec_line/advec_common_line.f90 b/HySoP/src/scalesInterface/particles/advec_line/advec_common_line.f90
index 7b4325354..f4757f4a4 100644
--- a/HySoP/src/scalesInterface/particles/advec_line/advec_common_line.f90
+++ b/HySoP/src/scalesInterface/particles/advec_line/advec_common_line.f90
@@ -87,8 +87,8 @@ subroutine AC_obtain_receivers_line(direction, ind_group, rece_ind_min, rece_ind
integer, dimension(2), intent(in) :: ind_group
integer, dimension(2), intent(out) :: rece_gap, send_gap
integer, dimension(MPI_STATUS_SIZE) :: statut
- ! Others
- integer :: proc_gap ! gap between a processus coordinate (along the current
+ ! Others
+ integer :: proc_gap ! gap between a processus coordinate (along the current
! direction) into the mpi-topology and my coordinate
integer :: rece_gapP ! gap between the coordinate of the previous processus (in the current direction)
! and the processes of maximal coordinate which will receive information from it
@@ -99,16 +99,16 @@ subroutine AC_obtain_receivers_line(direction, ind_group, rece_ind_min, rece_ind
integer :: send_request_bis ! mpi status of nonblocking send
integer :: ierr ! mpi error code
integer, dimension(2) :: tag_table ! some mpi message tag
- logical, dimension(:,:), allocatable:: test_request
- integer, dimension(:,:), allocatable:: s_request
+ logical, dimension(:,:), allocatable:: test_request
+ integer, dimension(:,:), allocatable:: s_request
tag_min = 5
tag_max = 6
send_gap = 3*mesh_sc%N(direction)
- rece_gap(1) = floor(real(rece_ind_min-1)/mesh_sc%N_proc(direction))
- rece_gap(2) = floor(real(rece_ind_max-1)/mesh_sc%N_proc(direction))
+ rece_gap(1) = floor(real(rece_ind_min-1, WP)/mesh_sc%N_proc(direction))
+ rece_gap(2) = floor(real(rece_ind_max-1, WP)/mesh_sc%N_proc(direction))
! ===== Communicate with my neigbors -> obtain ghost ! ====
! Compute their rank
@@ -130,7 +130,7 @@ subroutine AC_obtain_receivers_line(direction, ind_group, rece_ind_min, rece_ind
do proc_gap = rece_gap(1), rece_gap(2)
! Compute the rank of the target processus
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, rankN, ierr)
- ! Determine if I am the the first or the last processes (considering the current directory)
+ ! Determine if I am the the first or the last processes (considering the current directory)
! to require information from this processus
if (proc_gap>rece_gapP-1) then
if(rankN /= D_rank(direction)) then
@@ -180,9 +180,9 @@ end subroutine AC_obtain_receivers_line
!! @param[in,out] p_V = particle velocity (along the current direction)
!! @details
!! A RK2 scheme is used to advect the particles : the midlle point scheme. An
-!! intermediary position "p_pos_bis(i) = p_pos(i) + V(i)*dt/2" is computed and then
-!! the numerical velocity of each particles is computed as the interpolation of V in
-!! this point. This field is used to advect the particles at the seconde order in time :
+!! intermediary position "p_pos_bis(i) = p_pos(i) + V(i)*dt/2" is computed and then
+!! the numerical velocity of each particles is computed as the interpolation of V in
+!! this point. This field is used to advect the particles at the seconde order in time :
!! p_pos(t+dt, i) = p_pos(i) + p_V(i).
!! The group line indice is used to ensure using unicity of each mpi message tag.
subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
@@ -212,7 +212,7 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
integer :: rece_ind_max ! the maximal indice used in velocity interpolation
integer :: ind, ind_com ! indices
integer :: pos, pos_old ! indices of the mesh point wich preceed the particle position
- integer :: proc_gap, gap! distance between my (mpi) coordonate and coordinate of the
+ integer :: proc_gap, gap! distance between my (mpi) coordonate and coordinate of the
! processus associated to a given position
integer, dimension(:), allocatable :: rece_rank ! rank of processus wich send me information
integer :: send_rank ! rank of processus to wich I send information
@@ -277,7 +277,7 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank /= D_rank(direction)) then
! I - Receive messages about what I have to send
- ! Ia - Compute reception tag = concatenation of (rank+1), ind_group(1), ind_group(2), direction et unique Id.
+ ! Ia - Compute reception tag = concatenation of (rank+1), ind_group(1), ind_group(2), direction et unique Id.
tag = compute_tag(ind_group, tag_velo_range, direction, -proc_gap)
! Ib - Receive the message
call mpi_recv(send_range(1), 2, MPI_INTEGER, send_rank, tag, D_comm(direction), rece_status, ierr)
@@ -301,35 +301,35 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
! IIb - Receive message
gap = proc_gap*mesh_sc%N_proc(direction)
rece_range(1) = max(rece_ind_min, gap+1) ! fortran => indice start from 0
- rece_range(2) = min(rece_ind_max, gap+mesh_sc%N_proc(direction))
+ rece_range(2) = min(rece_ind_max, gap+mesh_sc%N_proc(direction))
msg_size = rece_range(2)-rece_range(1)+1
call mpi_Irecv(V_buffer(ind), msg_size, MPI_DOUBLE_PRECISION, rece_rank(proc_gap), tag, D_comm(direction), &
- & rece_request(proc_gap), ierr)
+ & rece_request(proc_gap), ierr)
ind = ind + msg_size
end if
end do
- ! -- Compute the interpolated velocity
+ ! -- Compute the interpolated velocity
! Compute the interpolation weight and update the pointers Vp and Vm
! Initialisation of reccurence process
ind = 1
pos = floor(p_pos_bis(ind))
weight(ind) = p_pos_bis(ind)-pos
! Vm = V(pos)
- proc_gap = floor(real(pos-1)/mesh_sc%N_proc(direction))
+ proc_gap = floor(real(pos-1, WP)/mesh_sc%N_proc(direction))
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank == D_rank(direction)) then
Vm(ind)%pter => p_V_bis(pos-proc_gap*mesh_sc%N_proc(direction))
- else
+ else
ind_com = ind_com + 1
Vm(ind)%pter => V_buffer(ind_com)
end if
! Vp = V(pos+1)
- proc_gap = floor(real(pos+1-1)/mesh_sc%N_proc(direction))
+ proc_gap = floor(real(pos+1-1, WP)/mesh_sc%N_proc(direction))
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank == D_rank(direction)) then
Vp(ind)%pter => p_V_bis(pos+1-proc_gap*mesh_sc%N_proc(direction))
- else
+ else
ind_com = ind_com + 1
Vp(ind)%pter => V_buffer(ind_com)
end if
@@ -348,31 +348,31 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
! The particle follows the previous one
Vm(ind)%pter => Vp(ind-1)%pter
! Vp = V(pos+1)
- proc_gap = floor(real(pos+1-1)/mesh_sc%N_proc(direction)) ! fortran -> indice starts from 1
+ proc_gap = floor(real(pos+1-1, WP)/mesh_sc%N_proc(direction)) ! fortran -> indice starts from 1
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank == D_rank(direction)) then
Vp(ind)%pter => p_V_bis(pos+1-proc_gap*mesh_sc%N_proc(direction))
- else
+ else
ind_com = ind_com + 1
Vp(ind)%pter => V_buffer(ind_com)
end if
case(2)
! pos = pos_old +2, wich correspond to "extention"
! Vm = V(pos)
- proc_gap = floor(real(pos-1)/mesh_sc%N_proc(direction))
+ proc_gap = floor(real(pos-1, WP)/mesh_sc%N_proc(direction))
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank == D_rank(direction)) then
Vm(ind)%pter => p_V_bis(pos-proc_gap*mesh_sc%N_proc(direction))
- else
+ else
ind_com = ind_com + 1
Vm(ind)%pter => V_buffer(ind_com)
end if
! Vp = V(pos+1)
- proc_gap = floor(real(pos+1-1)/mesh_sc%N_proc(direction))
+ proc_gap = floor(real(pos+1-1, WP)/mesh_sc%N_proc(direction))
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank == D_rank(direction)) then
Vp(ind)%pter => p_V_bis(pos+1-proc_gap*mesh_sc%N_proc(direction))
- else
+ else
ind_com = ind_com + 1
Vp(ind)%pter => V_buffer(ind_com)
end if
@@ -415,7 +415,7 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
end do
deallocate(s_request_bis)
- ! Deallocation
+ ! Deallocation
deallocate(rece_rank)
deallocate(rece_request)
deallocate(V_buffer)
@@ -434,7 +434,7 @@ end subroutine AC_particle_velocity_line
!! @param[in] ind_group = coordinate of the current group of lines
!! @param[in] p_V = particle velocity (along the current direction)
!! @param[out] bl_type = table of blocks type (center of left)
-!! @param[out] bl_tag = inform about tagged particles (bl_tag(ind_bl)=1 if the end of the bl_ind-th block
+!! @param[out] bl_tag = inform about tagged particles (bl_tag(ind_bl)=1 if the end of the bl_ind-th block
!! and the begining of the following one is tagged)
!! @details
!! This subroutine deals with a single line. For each line of this group, it
@@ -487,7 +487,7 @@ subroutine AC_type_and_block_line(dt, direction, ind_group, p_V, &
! Send message
call mpi_Isend(bl_lambdaMin(1), 1, MPI_DOUBLE_PRECISION, rankP, tag_table(1), D_comm(direction), send_request(1), ierr)
! Receive it
- call mpi_Irecv(lambN, 1, MPI_DOUBLE_PRECISION, rankN, tag_table(1), D_comm(direction), rece_request(1), ierr)
+ call mpi_Irecv(lambN, 1, MPI_DOUBLE_PRECISION, rankN, tag_table(1), D_comm(direction), rece_request(1), ierr)
! -- For the last block (1/2) --
! The processus contains only its first half => exchange ghost with the next processus
@@ -496,7 +496,7 @@ subroutine AC_type_and_block_line(dt, direction, ind_group, p_V, &
! Send message
call mpi_Isend(bl_lambdaMin(ind), 1, MPI_DOUBLE_PRECISION, rankN, tag_table(2), D_comm(direction), send_request(2), ierr)
! Receive it
- call mpi_Irecv(lambP, 1, MPI_DOUBLE_PRECISION, rankP, tag_table(2), D_comm(direction), rece_request(2), ierr)
+ call mpi_Irecv(lambP, 1, MPI_DOUBLE_PRECISION, rankP, tag_table(2), D_comm(direction), rece_request(2), ierr)
! -- For the "middle" block --
do ind = 2, bl_nb(direction)
@@ -554,10 +554,10 @@ end subroutine AC_type_and_block_line
!! Obtain the list of processus which contains some particles which belong to
!! my subdomains after their advection (and thus which will be remeshing into
!! my subdomain). This result is return as an interval [send_min; send_max].
-!! All the processus whose coordinate (into the current direction) belong to
+!! All the processus whose coordinate (into the current direction) belong to
!! this segment are involved into scalar remeshing into the current
!! subdomains. This routine does not involve any computation to determine if
-!! a processus is the first or the last processes (considering its coordinate along
+!! a processus is the first or the last processes (considering its coordinate along
!! the current directory) to send remeshing information to a given processes.
!! It directly compute it using contraints on velocity (as in corrected lambda
!! scheme) When possible use it rather than AC_obtain_senders_com
@@ -579,7 +579,7 @@ subroutine AC_obtain_senders_line(send_i_min, send_i_max, direction, ind_group,
integer, dimension(2), intent(out) :: rece_proc
integer, dimension(MPI_STATUS_SIZE) :: statut
! Other local variable
- integer(kind=4) :: proc_gap ! gap between a processus coordinate (along the current
+ integer(kind=4) :: proc_gap ! gap between a processus coordinate (along the current
! direction) into the mpi-topology and my coordinate
integer :: rankP, rankN ! processus rank for shift (P= previous, N = next)
integer, dimension(2) :: tag_table ! mpi message tag (for communicate rece_proc(1) and rece_proc(2))
@@ -590,8 +590,8 @@ subroutine AC_obtain_senders_line(send_i_min, send_i_max, direction, ind_group,
rece_proc = 3*mesh_sc%N(direction)
- proc_min = floor(real(send_i_min-1)/mesh_sc%N_proc(direction))
- proc_max = floor(real(send_i_max-1)/mesh_sc%N_proc(direction))
+ proc_min = floor(real(send_i_min-1, WP)/mesh_sc%N_proc(direction))
+ proc_max = floor(real(send_i_max-1, WP)/mesh_sc%N_proc(direction))
allocate(send_request(proc_min:proc_max,3))
send_request(:,3) = 0
@@ -667,7 +667,7 @@ end subroutine AC_obtain_senders_line
!! @param[in] send_buffer = buffer use to remesh the scalar before to send it to the right subdomain
!! @param[in,out] scal1D = mono-dimensionnal scalar field to advect
!! @details
-!! Remeshing are done in a local buffer. This subroutine distribute this buffer
+!! Remeshing are done in a local buffer. This subroutine distribute this buffer
!! to the right processes, receive the buffer send and update the scalar field.
subroutine AC_bufferToScalar_line(direction, ind_group, send_i_min, send_i_max, proc_min, proc_max, rece_proc,send_buffer, scal1D)
@@ -680,7 +680,7 @@ subroutine AC_bufferToScalar_line(direction, ind_group, send_i_min, send_i_max,
integer, dimension(2), intent(in) :: ind_group
integer, intent(in) :: send_i_min
integer, intent(in) :: send_i_max
- integer, dimension(2), intent(in) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the
+ integer, dimension(2), intent(in) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the
! one from which I will receive data
integer, intent(in) :: proc_min ! smaller gap between me and the processes to where I send data
integer, intent(in) :: proc_max ! smaller gap between me and the processes to where I send data
@@ -690,7 +690,7 @@ subroutine AC_bufferToScalar_line(direction, ind_group, send_i_min, send_i_max,
! Variables used to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
! design something I send (resp. I receive).
integer :: i ! table indice
- integer :: proc_gap ! gap between my Y-coordinate and the one of the processus
+ integer :: proc_gap ! gap between my Y-coordinate and the one of the processus
real(WP), dimension(:), allocatable :: rece_buffer ! buffer use to stock received scalar field
integer :: send_gap ! number of mesh between my and another processus
integer,dimension(:,:), allocatable :: rece_range ! range of (local) indice where the received scalar field has to be save
@@ -709,7 +709,7 @@ subroutine AC_bufferToScalar_line(direction, ind_group, send_i_min, send_i_max,
integer :: tag ! mpi message tag
! with wich I communicate.
- ! ===== Receive information =====
+ ! ===== Receive information =====
! -- Allocate field --
allocate(rece_rank(rece_proc(1):rece_proc(2)))
allocate(rece_range(2,rece_proc(1):rece_proc(2))) ! be careful that mpi use contiguous memory element
@@ -733,7 +733,7 @@ subroutine AC_bufferToScalar_line(direction, ind_group, send_i_min, send_i_max,
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
send_gap = proc_gap*mesh_sc%N_proc(direction)
send_range(1, proc_gap) = max(send_i_min, send_gap+1) ! fortran => indice start from 0
- send_range(2, proc_gap) = min(send_i_max, send_gap+mesh_sc%N_proc(direction))
+ send_range(2, proc_gap) = min(send_i_max, send_gap+mesh_sc%N_proc(direction))
if (send_rank/=D_rank(direction)) then
! Determine quantity of information to send
comm_size = send_range(2, proc_gap)-send_range(1, proc_gap)+1
diff --git a/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_formula_line.f90 b/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_formula_line.f90
index 2cca2e257..8bd63c646 100644
--- a/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_formula_line.f90
+++ b/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_formula_line.f90
@@ -6,14 +6,14 @@
! MODULE: advec_remeshing_line
!
!
-! DESCRIPTION:
+! DESCRIPTION:
!> This module gathers all the remeshing formula. These interpolation
!!polynom allow to re-distribute particles on mesh grid at each
!! iterations. - old version for advection without group of line
!! @details
!! It provides lambda 2 corrected, lambda 4 corrected and M'6 remeshing formula.
!! The remeshing of type "lambda corrected" are design for large time
-!! step. The M'6 formula appears as being stable for large time step, but
+!! step. The M'6 formula appears as being stable for large time step, but
!! the numerical analysis remains todo.
!! This module also provide some wraper to remesh a complete line
!! of particles (with the different formula) and to do it either on a
@@ -47,7 +47,7 @@ contains
!! @param[in] ind_max = maximal indice of the send buffer
!! @param[in, out] send_buffer = buffer use to remesh the scalar before to send it to the right subdomain
!! @details
-!! Use corrected lambda 2 remeshing formula.
+!! Use corrected lambda 2 remeshing formula.
!! This remeshing formula depends on the particle type :
!! 1 - Is the particle tagged ?
!! 2 - Does it belong to a centered or a left block ?
@@ -126,7 +126,7 @@ end subroutine AC_remesh_lambda2corrected_basic
!! @param[in] ind_max = maximal indice of the send buffer
!! @param[in, out] send_buffer = buffer use to remesh the scalar before to send it to the right subdomain
!! @details
-!! Use corrected lambda 2 remeshing formula.
+!! Use corrected lambda 2 remeshing formula.
!! This remeshing formula depends on the particle type :
!! 1 - Is the particle tagged ?
!! 2 - Does it belong to a centered or a left block ?
@@ -157,7 +157,7 @@ subroutine AC_remesh_lambda4corrected_basic(direction, p_pos_adim, scal1D, bl_ty
do p_ind = 1, mesh_sc%N_proc(direction), bl_size
bl_ind = p_ind/bl_size + 1
- if (bl_tag(bl_ind)) then
+ if (bl_tag(bl_ind)) then
! Tagged case
if (bl_type(bl_ind)) then
! tagged, the first particle belong to a centered block and the last to left block.
@@ -170,14 +170,14 @@ subroutine AC_remesh_lambda4corrected_basic(direction, p_pos_adim, scal1D, bl_ty
end if
else
! No tag
- if (bl_type(bl_ind)) then
+ if (bl_type(bl_ind)) then
call AC_remesh_O4_center(p_pos_adim(p_ind),scal1D(p_ind), send_buffer)
call AC_remesh_O4_center(p_pos_adim(p_ind+1),scal1D(p_ind+1), send_buffer)
else
call AC_remesh_O4_left(p_pos_adim(p_ind),scal1D(p_ind), send_buffer)
call AC_remesh_O4_left(p_pos_adim(p_ind+1),scal1D(p_ind+1), send_buffer)
end if
- if (bl_type(bl_ind+1)) then
+ if (bl_type(bl_ind+1)) then
call AC_remesh_O4_center(p_pos_adim(p_ind+2),scal1D(p_ind+2), send_buffer)
call AC_remesh_O4_center(p_pos_adim(p_ind+3),scal1D(p_ind+3), send_buffer)
else
@@ -198,7 +198,7 @@ subroutine AC_remesh_left(pos_adim, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
-
+
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
real(WP), dimension(send_j_min:send_j_max), intent(inout) :: buffer
@@ -209,9 +209,9 @@ subroutine AC_remesh_left(pos_adim, sca, buffer)
! Mesh point used in remeshing formula
j0 = floor(pos_adim)
!j0 = floor(pos/d_sc(2))
-
+
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
bM=0.5*y0*(y0-1.)
@@ -248,7 +248,7 @@ subroutine AC_remesh_center(pos_adim, sca, buffer)
!j0 = nint(pos/d_sc(2))
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
!y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
! Interpolation weights
@@ -261,7 +261,7 @@ subroutine AC_remesh_center(pos_adim, sca, buffer)
buffer(j0-1) = buffer(j0-1) + bM*sca
buffer(j0) = buffer(j0) + b0*sca
buffer(j0+1) = buffer(j0+1) + bP*sca
-
+
end subroutine AC_remesh_center
@@ -271,9 +271,9 @@ end subroutine AC_remesh_center
!! @param[in] posP_ad = adimensionned position of the second particle
!! @param[in] scaP = scalar advected by this particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-!! @details
+!! @details
!! Remeshing formula devoted to tagged particles.
-!! The particle group send into argument is composed of a block end and of the
+!! The particle group send into argument is composed of a block end and of the
!! begining of the next block. The first particles belong to a centered block
!! and the last to a left one. The block have difference indice (tagged
!! particles) and we have to use corrected formula.
@@ -299,9 +299,9 @@ subroutine AC_remesh_tag_CL(pos_adim, sca, posP_ad, scaP, buffer)
jM=j0-1
jP=j0+1
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
!y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
- y0_bis = (posP_ad - dble(j0_bis))
+ y0_bis = (posP_ad - real(j0_bis, WP))
!y0_bis = (posP - dble(j0_bis)*d_sc(2))/d_sc(2)
aM=0.5*y0*(y0-1)
@@ -323,9 +323,9 @@ end subroutine AC_remesh_tag_CL
!! @param[in] posP_ad = adimensionned position of the second particle
!! @param[in] scaP = scalar advected by this particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-!! @details
+!! @details
!! Remeshing formula devoted to tagged particles.
-!! The particle group send into argument is composed of a block end and of the
+!! The particle group send into argument is composed of a block end and of the
!! begining of the next block. The first particles belong to a left block
!! and the last to a centered one. The block have difference indice (tagged
!! particles) and we have to use corrected formula.
@@ -354,13 +354,13 @@ subroutine AC_remesh_tag_LC(pos_adim, sca, posP_ad, scaP, buffer)
jP=j0+1
jP2=j0+2
jP3=j0+3
-
+
! Distance to mesh point
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
!y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
- y0_bis = (posP_ad - dble(j0_bis))
+ y0_bis = (posP_ad - real(j0_bis, WP))
!y0_bis = (posP - dble(j0_bis)*d_sc(2))/d_sc(2)
-
+
! Interpolation weight
a0=1-y0**2
aP=y0
@@ -391,7 +391,7 @@ subroutine AC_remesh_O4_left(pos_adim, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
-
+
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
real(WP), dimension(send_j_min:send_j_max), intent(inout) :: buffer
@@ -402,9 +402,9 @@ subroutine AC_remesh_O4_left(pos_adim, sca, buffer)
! Mesh point used in remeshing formula
j0 = floor(pos_adim)
!j0 = floor(pos/d_sc(2))
-
+
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
!bM2=(y0-2.)*(y0-1.)*y0*(y0+1.)/24.0
@@ -446,9 +446,9 @@ subroutine AC_remesh_O4_center(pos_adim, sca, buffer)
real(WP) :: y0 ! adimensionned distance to mesh points
! Mesh point used in remeshing formula
j0 = nint(pos_adim)
-
+
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
!bM2=(y0-2.)*(y0-1.)*y0*(y0+1.)/24.0
@@ -468,7 +468,7 @@ subroutine AC_remesh_O4_center(pos_adim, sca, buffer)
buffer(j0) = buffer(j0) + b0*sca
buffer(j0+1) = buffer(j0+1) + bP*sca
buffer(j0+2) = buffer(j0+2) + bP2*sca
-
+
end subroutine AC_remesh_O4_center
@@ -483,9 +483,9 @@ end subroutine AC_remesh_O4_center
!! @param[in] posP2_ad= adimensionned position of the fourth (and last) particle
!! @param[in] scaP2 = scalar advected by this particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-!! @details
+!! @details
!! Remeshing formula devoted to tagged particles.
-!! The particle group send into argument is composed of a block end and of the
+!! The particle group send into argument is composed of a block end and of the
!! begining of the next block. The first particles belong to a centered block
!! and the last to a left one. The block have difference indice (tagged
!! particles) and we have to use corrected formula.
@@ -514,10 +514,10 @@ subroutine AC_remesh_O4_tag_CL(posM_ad, scaM, pos_adim, sca, posP_ad, scaP, posP
jP2= floor(posP2_ad)
! Distance to mesh point
- yM = (posM_ad - dble(jM))
- y0 = (pos_adim - dble(j0))
- yP = (posP_ad - dble(jP))
- yP2= (posP2_ad - dble(jP2))
+ yM = (posM_ad - real(jM, WP))
+ y0 = (pos_adim - real(j0, WP))
+ yP = (posP_ad - real(jP, WP))
+ yP2= (posP2_ad - real(jP2, WP))
! Interpolation weights
!aM3=(yM-2.)*(yM-1.)*yM*(yM+1.)/24.0
@@ -533,8 +533,8 @@ subroutine AC_remesh_O4_tag_CL(posM_ad, scaM, pos_adim, sca, posP_ad, scaP, posP
bM2=y0*(2.+y0*(-1.+y0*(-2.+y0)))/24.0
!bM =(2.-y0)*(y0-1.)*y0*(y0+2.)/6.0
bM =y0*(-4.+y0*(4.+y0*(1.-y0)))/6.0
- !bP =((y0+1)-1.)*(y0+1)*((y0+1)+1.)*((y0+1)+2.)/24.0
- bP =y0*(6.+y0*(11+y0*(6+y0)))/24.0
+ !bP =((y0+1)-1.)*(y0+1)*((y0+1)+1.)*((y0+1)+2.)/24.0
+ bP =y0*(6.+y0*(11+y0*(6+y0)))/24.0
!b0 =((y0-2.)*(y0-1.)*(y0+1.)*(y0+2.)/4.0) + ((2.-y0)*y0*(y0+1.)*(y0+2.)/6.0) &
! & + ((y0-1.)*y0*(y0+1.)*(y0+2.)/24.0) - bP
b0 = 1. - (bM2+bM+bP)
@@ -559,7 +559,7 @@ subroutine AC_remesh_O4_tag_CL(posM_ad, scaM, pos_adim, sca, posP_ad, scaP, posP
e0 = 1. - (eP+eP2+eP3)
! remeshing
- buffer(j0-3) = buffer(j0-3) +aM3*scaM
+ buffer(j0-3) = buffer(j0-3) +aM3*scaM
buffer(j0-2) = buffer(j0-2) +aM2*scaM +bM2*sca
buffer(j0-1) = buffer(j0-1) + aM*scaM + bM*sca + cM*scaP
buffer(j0) = buffer(j0) + a0*scaM + b0*sca + c0*scaP + e0*scaP2
@@ -582,9 +582,9 @@ end subroutine AC_remesh_O4_tag_CL
!! @param[in] posP2_ad= adimensionned position of the fourth (and last) particle
!! @param[in] scaP2 = scalar advected by this particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-!! @details
+!! @details
!! Remeshing formula devoted to tagged particles.
-!! The particle group send into argument is composed of a block end and of the
+!! The particle group send into argument is composed of a block end and of the
!! begining of the next block. The first particles belong to a left block
!! and the last to a centered one. The block have difference indice (tagged
!! particles) and we have to use corrected formula.
@@ -614,10 +614,10 @@ subroutine AC_remesh_O4_tag_LC(posM_ad, scaM, pos_adim, sca, posP_ad, scaP, posP
jP2= nint(posP2_ad)
! Distance to mesh point
- yM = (posM_ad - dble(jM))
- y0 = (pos_adim - dble(j0))
- yP = (posP_ad - dble(jP))
- yP2= (posP2_ad - dble(jP2))
+ yM = (posM_ad - real(jM, WP))
+ y0 = (pos_adim - real(j0, WP))
+ yP = (posP_ad - real(jP, WP))
+ yP2= (posP2_ad - real(j0, WP))
! Interpolation weights
!aM3=(yM-2.)*(yM-1.)*yM*(yM+1.)/24.0
@@ -626,7 +626,7 @@ subroutine AC_remesh_O4_tag_LC(posM_ad, scaM, pos_adim, sca, posP_ad, scaP, posP
aM2 =yM*(-4.+yM*(4.+yM*(1.-yM)))/6.0
!aM =(yM-2.)*(yM-1.)*(yM+1.)*(yM+2.)/4.0
aM =(4.+(yM**2)*(-5.+yM**2))/4.0
- !a0 =((2.-yM)*yM*(yM+1.)*(yM+2.)/6.0)
+ !a0 =((2.-yM)*yM*(yM+1.)*(yM+2.)/6.0)
a0 =yM*(4+yM*(4-yM*(1.+yM)))/6.0
!aP2=(((yM-1.)-1.)*(yM-1.)*((yM-1.)+1.)*((yM-1.)+2.)/24.0)
!aP2=yM*(yM-2.)*(yM-1.)*(yM+1.)/24.0
@@ -683,7 +683,7 @@ subroutine AC_remesh_O4_tag_LC(posM_ad, scaM, pos_adim, sca, posP_ad, scaP, posP
eP = 1.0 - (e0+eP2+eP3+eP4+eP5)
! remeshing
- buffer(j0-3) = buffer(j0-3) +aM3*scaM
+ buffer(j0-3) = buffer(j0-3) +aM3*scaM
buffer(j0-2) = buffer(j0-2) +aM2*scaM +bM2*sca
buffer(j0-1) = buffer(j0-1) + aM*scaM + bM*sca + cM*scaP
buffer(j0) = buffer(j0) + a0*scaM + b0*sca + c0*scaP + e0*scaP2
@@ -704,7 +704,7 @@ subroutine AC_remesh_Mprime6(pos_adim, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
-
+
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
real(WP), dimension(send_j_min:send_j_max), intent(inout) :: buffer
@@ -715,9 +715,9 @@ subroutine AC_remesh_Mprime6(pos_adim, sca, buffer)
! Mesh point used in remeshing formula
j0 = floor(pos_adim)
-
+
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
!bM2 =-(((y0+2.)-2)*(5.*(y0+2.)-8.)*((y0+2.)-3.)**3)/24.
diff --git a/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_line.F90 b/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_line.F90
index 6331ae32d..68fdc2feb 100644
--- a/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_line.F90
+++ b/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_line.F90
@@ -187,7 +187,7 @@ subroutine Xremesh_O4(direction, ind_group, gs, p_pos_adim, p_V, j,k, scal, dt)
! ... and to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
! designes something I send (resp. I receive).
real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
- integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
! I will receive data
integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
@@ -456,7 +456,7 @@ subroutine Yremesh_O4(direction, ind_group, gs, p_pos_adim, p_V, i,k,scal, dt)
! ... and to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
! designes something I send (resp. I receive).
real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
- integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
! I will receive data
integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
@@ -723,7 +723,7 @@ subroutine Zremesh_O4(direction, ind_group, gs, p_pos_adim, p_V,i,j,scal, dt)
! ... and to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
! designes something I send (resp. I receive).
real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
- integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
! I will receive data
integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
@@ -939,8 +939,8 @@ subroutine AC_obtain_senders_group(direction, gp_s, ind_group, send_min, send_ma
rece_proc = 3*mesh_sc%N(direction)
max_size = size(rece_buffer) + 1
- proc_min = floor(real(send_min-1)/mesh_sc%N_proc(direction))
- proc_max = floor(real(send_max-1)/mesh_sc%N_proc(direction))
+ proc_min = floor(real(send_min-1, WP)/mesh_sc%N_proc(direction))
+ proc_max = floor(real(send_max-1, WP)/mesh_sc%N_proc(direction))
proc_min_abs = minval(proc_min)
proc_max_abs = maxval(proc_max)
@@ -1125,7 +1125,7 @@ subroutine AC_obtain_senders_com(direction, gp_s, ind_group, send_min, send_max,
integer, dimension(:,:), intent(in) :: send_min ! distance between me and processus wich send me information
integer, dimension(:,:), intent(in) :: send_max ! distance between me and processus wich send me information
! Other local variable
- integer(kind=4) :: proc_gap ! gap between a processus coordinate (along the current
+ integer(kind=4) :: proc_gap ! gap between a processus coordinate (along the current
! direction) into the mpi-topology and my coordinate
integer :: rankP, rankN ! processus rank for shift (P= previous, N = next)
integer, dimension(2) :: tag_table ! mpi message tag (for communicate rece_proc(1) and rece_proc(2))
@@ -1148,8 +1148,8 @@ subroutine AC_obtain_senders_com(direction, gp_s, ind_group, send_min, send_max,
rece_proc = 3*mesh_sc%N(direction)
max_size = size(rece_buffer) + 1
- proc_min = floor(real(send_min-1)/mesh_sc%N_proc(direction))
- proc_max = floor(real(send_max-1)/mesh_sc%N_proc(direction))
+ proc_min = floor(real(send_min-1, WP)/mesh_sc%N_proc(direction))
+ proc_max = floor(real(send_max-1, WP)/mesh_sc%N_proc(direction))
proc_min_abs = minval(proc_min)
proc_max_abs = maxval(proc_max)
diff --git a/HySoP/src/scalesInterface/particles/advec_remesh_Mprime.f90 b/HySoP/src/scalesInterface/particles/advec_remesh_Mprime.f90
index 2a18effc6..9725e81b2 100644
--- a/HySoP/src/scalesInterface/particles/advec_remesh_Mprime.f90
+++ b/HySoP/src/scalesInterface/particles/advec_remesh_Mprime.f90
@@ -176,7 +176,7 @@ subroutine AC_remesh_Mprime4_array(dir, pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
!bM = ((2.-(y0+1.))**2 * (1.-(y0+1.)))/2.
@@ -204,9 +204,10 @@ end subroutine AC_remesh_Mprime4_array
!> M'4 remeshing formula. - version for array of pointer
!! @author Chloe Mimeau, LJK
!! @param[in] pos_adim= adimensionned particle position
+!! @param[in] translat= translation to convert adimensionned particle position to the proper array index
!! @param[in] sca = scalar advected by the particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-subroutine AC_remesh_Mprime4_pter(pos_adim, sca, buffer)
+subroutine AC_remesh_Mprime4_pter(pos_adim, translat, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
@@ -214,6 +215,7 @@ subroutine AC_remesh_Mprime4_pter(pos_adim, sca, buffer)
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
type(real_pter), dimension(:), intent(inout) :: buffer
+ integer, intent(in) :: translat
! Other local variables
integer :: j0 ! indice of the nearest mesh points
real(WP) :: bM, b0, bP, bP2 ! interpolation weight for the particles
@@ -223,7 +225,10 @@ subroutine AC_remesh_Mprime4_pter(pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
+
+ ! translation to obtain the array index
+ j0 = j0 + translat
! Interpolation weights
!bM = ((2.-(y0+1.))**2 * (1.-(y0+1.)))/2.
@@ -270,7 +275,7 @@ subroutine AC_remesh_Mprime4_diff_array(dir, pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Compute coefficient for diffusion part
diff1 = 1.5*(1.-0.5*4.*sc_diff_dt_dx(sc_remesh_ind,dir))
@@ -305,9 +310,10 @@ end subroutine AC_remesh_Mprime4_diff_array
!! @author Jean-baptiste Lagaert, LEGI
!! @param[in] diff = diffusivity
!! @param[in] pos_adim= adimensionned particle position
+!! @param[in] translat= translation to convert adimensionned particle position to the proper array index
!! @param[in] sca = scalar advected by the particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-subroutine AC_remesh_Mprime4_diff_pter(pos_adim, sca, buffer)
+subroutine AC_remesh_Mprime4_diff_pter(pos_adim, translat, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
@@ -315,6 +321,7 @@ subroutine AC_remesh_Mprime4_diff_pter(pos_adim, sca, buffer)
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
type(real_pter), dimension(:), intent(inout) :: buffer
+ integer, intent(in) :: translat
! Other local variables
integer :: j0 ! indice of the nearest mesh points
real(WP) :: bM, b0, bP, bP2 ! interpolation weight for the particles
@@ -325,7 +332,10 @@ subroutine AC_remesh_Mprime4_diff_pter(pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
+
+ ! translation to obtain the array index
+ j0 = j0 + translat
! Compute coefficient for diffusion part
diff1 = 1.5*(1.-0.5*4.*sc_diff_dt_dx(sc_remesh_ind,current_dir))
@@ -375,7 +385,7 @@ subroutine AC_remesh_Mstar6_array(dir, pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
!bM2 =-(((y0+2.)-2)*(5.*(y0+2.)-8.)*((y0+2.)-3.)**3)/24.
@@ -413,9 +423,10 @@ end subroutine AC_remesh_Mstar6_array
!! determining order). - version for array of pointer
!! @author Jean-Baptiste Lagaert, LEGI/LJK
!! @param[in] pos_adim= adimensionned particle position
+!! @param[in] translat= translation to convert adimensionned particle position to the proper array index
!! @param[in] sca = scalar advected by the particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-subroutine AC_remesh_Mstar6_pter(pos_adim, sca, buffer)
+subroutine AC_remesh_Mstar6_pter(pos_adim, translat, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
@@ -423,6 +434,7 @@ subroutine AC_remesh_Mstar6_pter(pos_adim, sca, buffer)
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
type(real_pter), dimension(:), intent(inout) :: buffer
+ integer, intent(in) :: translat
! Other local variables
integer :: j0 ! indice of the nearest mesh points
real(WP) :: bM, bM2, b0, bP, bP2, bP3! interpolation weight for the particles
@@ -432,7 +444,10 @@ subroutine AC_remesh_Mstar6_pter(pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
+
+ ! translation to obtain the array index
+ j0 = j0 + translat
! Interpolation weights
!bM2 =-(((y0+2.)-2)*(5.*(y0+2.)-8.)*((y0+2.)-3.)**3)/24.
@@ -449,6 +464,8 @@ subroutine AC_remesh_Mstar6_pter(pos_adim, sca, buffer)
!b0 = (12. + y0**2*(-15. + y0*(-35. + (63. - 25.*y0)*y0)))/12.
b0 = 1. - (bM2+bM+bP+bP2+bP3)
+ !print *, j0, pos_adim
+
! remeshing
buffer(j0-2)%pter = buffer(j0-2)%pter + sca*bM2
buffer(j0-1)%pter = buffer(j0-1)%pter + sca*bM
@@ -484,7 +501,7 @@ subroutine AC_remesh_L4_4_array(dir, pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
bM2 = (y0*(y0*(y0*(y0*(y0*(y0*(y0*(y0*(-46. * y0 + 207.) - 354.) + 273.) - 80.) + 1.) - 2.)- 1.) + 2.)) / 24.
@@ -515,9 +532,10 @@ end subroutine AC_remesh_L4_4_array
!> Lambda(4,4) uncorrected remeshing formula (order 4 everywhere) - version for array of pointer
!! @author Jean-Baptiste Lagaert, LEGI/LJK
!! @param[in] pos_adim= adimensionned particle position
+!! @param[in] translat= translation to convert adimensionned particle position to the proper array index
!! @param[in] sca = scalar advected by the particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-subroutine AC_remesh_L4_4_pter(pos_adim, sca, buffer)
+subroutine AC_remesh_L4_4_pter(pos_adim, translat, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
@@ -525,6 +543,7 @@ subroutine AC_remesh_L4_4_pter(pos_adim, sca, buffer)
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
type(real_pter), dimension(:), intent(inout) :: buffer
+ integer, intent(in) :: translat
! Other local variables
integer :: j0 ! indice of the nearest mesh points
real(WP) :: bM, bM2, b0, bP, bP2, bP3! interpolation weight for the particles
@@ -534,7 +553,10 @@ subroutine AC_remesh_L4_4_pter(pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
+
+ ! translation to obtain the array index
+ j0 = j0 + translat
! Interpolation weights
bM2 = (y0*(y0*(y0*(y0*(y0*(y0*(y0*(y0*(-46. * y0 + 207.) - 354.) + 273.) - 80.) + 1.) - 2.)- 1.) + 2.)) / 24.
@@ -582,7 +604,7 @@ subroutine AC_remesh_L6_4_array(dir, pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
bM3 = (y0*(y0*(y0*(y0*(y0*(y0*(y0*(y0*(290. * y0 - 1305.) + 2231.) &
@@ -628,9 +650,10 @@ end subroutine AC_remesh_L6_4_array
!! - version for array of pointer
!! @author Chloe Mimeau, LJK
!! @param[in] pos_adim= adimensionned particle position
+!! @param[in] translat= translation to convert adimensionned particle position to the proper array index
!! @param[in] sca = scalar advected by the particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-subroutine AC_remesh_L6_4_pter(pos_adim, sca, buffer)
+subroutine AC_remesh_L6_4_pter(pos_adim, translat, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
@@ -638,6 +661,7 @@ subroutine AC_remesh_L6_4_pter(pos_adim, sca, buffer)
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
type(real_pter), dimension(:), intent(inout) :: buffer
+ integer, intent(in) :: translat
! Other local variables
integer :: j0 ! indice of the nearest mesh points
real(WP) :: bM3, bM2, bM, b0 ! interpolation weight for the particles
@@ -648,7 +672,10 @@ subroutine AC_remesh_L6_4_pter(pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
+
+ ! translation to obtain the array index
+ j0 = j0 + translat
! Interpolation weights
bM3 = (y0*(y0*(y0*(y0*(y0*(y0*(y0*(y0*(290. * y0 - 1305.) + 2231.) &
@@ -707,7 +734,7 @@ subroutine AC_remesh_L6_6_array(dir, pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
bM3 = (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (3604. * y0 - 23426.) + 63866.) &
@@ -759,9 +786,10 @@ end subroutine AC_remesh_L6_6_array
!> Lambda(6,6) uncorrected remeshing formula (order 6 everywhere) - version for array of pointer
!! @author Jean-Baptiste Lagaert, LEGI/LJK
!! @param[in] pos_adim= adimensionned particle position
+!! @param[in] translat= translation to convert adimensionned particle position to the proper array index
!! @param[in] sca = scalar advected by the particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-subroutine AC_remesh_L6_6_pter(pos_adim, sca, buffer)
+subroutine AC_remesh_L6_6_pter(pos_adim, translat, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
@@ -769,6 +797,7 @@ subroutine AC_remesh_L6_6_pter(pos_adim, sca, buffer)
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
type(real_pter), dimension(:), intent(inout) :: buffer
+ integer, intent(in) :: translat
! Other local variables
integer :: j0 ! indice of the nearest mesh points
real(WP) :: bM, bM2, bM3, b0, bP, bP2, bP3, bP4! interpolation weight for the particles
@@ -778,7 +807,10 @@ subroutine AC_remesh_L6_6_pter(pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
+
+ ! translation to obtain the array index
+ j0 = j0 + translat
! Interpolation weights
bM3 = (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (y0 * (3604. * y0 - 23426.) + 63866.) &
@@ -845,7 +877,7 @@ subroutine AC_remesh_L8_4_array(dir, pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
bM4 = (y0*(y0*(y0*(y0*(y0*(y0*(y0*(y0*(-3569. * y0 + 16061.) &
@@ -908,9 +940,10 @@ end subroutine AC_remesh_L8_4_array
!! - version for array of pointer
!! @author Chloe Mimeau, LJK
!! @param[in] pos_adim= adimensionned particle position
+!! @param[in] translat= translation to convert adimensionned particle position to the proper array index
!! @param[in] sca = scalar advected by the particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-subroutine AC_remesh_L8_4_pter(pos_adim, sca, buffer)
+subroutine AC_remesh_L8_4_pter(pos_adim, translat, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
@@ -918,6 +951,7 @@ subroutine AC_remesh_L8_4_pter(pos_adim, sca, buffer)
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
type(real_pter), dimension(:), intent(inout) :: buffer
+ integer, intent(in) :: translat
! Other local variables
integer :: j0 ! indice of the nearest mesh points
real(WP) :: bM4, bM3, bM2, bM, b0 ! interpolation weight for the particles
@@ -928,7 +962,10 @@ subroutine AC_remesh_L8_4_pter(pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
+
+ ! translation to obtain the array index
+ j0 = j0 + translat
! Interpolation weights
bM4 = (y0*(y0*(y0*(y0*(y0*(y0*(y0*(y0*(-3569. * y0 + 16061.) &
@@ -1000,7 +1037,7 @@ subroutine AC_remesh_Mprime8_array(dir, pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
! M'8 = 15/8*M8 + 9/8 * x * M'8 + 1/8 * x^2 * M''8
@@ -1069,9 +1106,10 @@ end subroutine AC_remesh_Mprime8_array
!! @author Jean-Baptiste Lagaert, LEGI/LJK
!! @param[in] dir = current direction
!! @param[in] pos_adim= adimensionned particle position
+!! @param[in] translat= translation to convert adimensionned particle position to the proper array index
!! @param[in] sca = scalar advected by the particle
!! @param[in,out] buffer = temporaly remeshed scalar field
-subroutine AC_remesh_Mprime8_pter(pos_adim, sca, buffer)
+subroutine AC_remesh_Mprime8_pter(pos_adim, translat, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
@@ -1079,6 +1117,7 @@ subroutine AC_remesh_Mprime8_pter(pos_adim, sca, buffer)
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
type(real_pter), dimension(:), intent(inout) :: buffer
+ integer, intent(in) :: translat
! Other local variables
integer :: j0 ! indice of the nearest mesh points
real(WP) :: y0 ! adimensionned distance to mesh points
@@ -1088,7 +1127,10 @@ subroutine AC_remesh_Mprime8_pter(pos_adim, sca, buffer)
j0 = floor(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
+
+ ! translation to obtain the array index
+ j0 = j0 + translat
! Interpolation weights
! M'8 = 15/8*M8 + 9/8 * x * M'8 + 1/8 * x^2 * M''8
diff --git a/HySoP/src/scalesInterface/particles/advec_remesh_lambda.f90 b/HySoP/src/scalesInterface/particles/advec_remesh_lambda.f90
index 400375f38..f93168dcd 100644
--- a/HySoP/src/scalesInterface/particles/advec_remesh_lambda.f90
+++ b/HySoP/src/scalesInterface/particles/advec_remesh_lambda.f90
@@ -679,7 +679,7 @@ subroutine AC_remesh_O2_array(dir, pos_adim, sca, bl_type, buffer)
end if
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
bM=0.5*y0*(y0-1.)
@@ -727,7 +727,7 @@ subroutine AC_remesh_O2_pter(pos_adim, sca, bl_type, buffer)
end if
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
bM=0.5*y0*(y0-1.)
@@ -777,10 +777,10 @@ subroutine AC_remesh_tag_CL_array(dir, pos_adim, sca, posP_ad, scaP, buffer)
j0_bis = floor(posP_ad)
!j0_bis = floor(posP/d_sc(2))
- y0 = (pos_adim - dble(j0))
- !y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
- y0_bis = (posP_ad - dble(j0_bis))
- !y0_bis = (posP - dble(j0_bis)*d_sc(2))/d_sc(2)
+ y0 = (pos_adim - real(j0, WP))
+ !y0 = (pos - real(j0, WP)*d_sc(2))/d_sc(2)
+ y0_bis = (posP_ad - real(j0_bis, WP))
+ !y0_bis = (posP - real(j0_bis, WP)*d_sc(2))/d_sc(2)
aM=0.5*y0*(y0-1)
a0=1.-aM
@@ -832,10 +832,10 @@ subroutine AC_remesh_tag_CL_pter(pos_adim, sca, posP_ad, scaP, buffer)
jM=j0-1
jP=j0+1
- y0 = (pos_adim - dble(j0))
- !y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
- y0_bis = (posP_ad - dble(j0_bis))
- !y0_bis = (posP - dble(j0_bis)*d_sc(2))/d_sc(2)
+ y0 = (pos_adim - real(j0, WP))
+ !y0 = (pos - real(j0, WP)*d_sc(2))/d_sc(2)
+ y0_bis = (posP_ad - real(j0_bis, WP))
+ !y0_bis = (posP - real(j0_bis, WP)*d_sc(2))/d_sc(2)
aM=0.5*y0*(y0-1)
a0=1.-aM
@@ -889,13 +889,13 @@ subroutine AC_remesh_tag_LC_array(dir, pos_adim, sca, posP_ad, scaP, buffer)
jP=j0+1
jP2=j0+2
jP3=j0+3
-
+
! Distance to mesh point
- y0 = (pos_adim - dble(j0))
- !y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
- y0_bis = (posP_ad - dble(j0_bis))
- !y0_bis = (posP - dble(j0_bis)*d_sc(2))/d_sc(2)
-
+ y0 = (pos_adim - real(j0, WP))
+ !y0 = (pos - real(j0, WP)*d_sc(2))/d_sc(2)
+ y0_bis = (posP_ad - real(j0_bis, WP))
+ !y0_bis = (posP - real(j0_bis, WP)*d_sc(2))/d_sc(2)
+
! Interpolation weight
a0=1-y0**2
aP=y0
@@ -960,13 +960,13 @@ subroutine AC_remesh_tag_LC_pter(pos_adim, sca, posP_ad, scaP, buffer)
jP=j0+1
jP2=j0+2
jP3=j0+3
-
+
! Distance to mesh point
- y0 = (pos_adim - dble(j0))
- !y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
- y0_bis = (posP_ad - dble(j0_bis))
- !y0_bis = (posP - dble(j0_bis)*d_sc(2))/d_sc(2)
-
+ y0 = (pos_adim - real(j0, WP))
+ !y0 = (pos - real(j0, WP)*d_sc(2))/d_sc(2)
+ y0_bis = (posP_ad - real(j0_bis, WP))
+ !y0_bis = (posP - real(j0_bis, WP)*d_sc(2))/d_sc(2)
+
! Interpolation weight
a0=1-y0**2
aP=y0
@@ -1032,7 +1032,7 @@ subroutine AC_remesh_limitO2_array(dir, pos_adim, sca, bl_type, limit, buffer)
end if
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
bM=0.5*((y0-0.5)**2) - limit(1)
@@ -1090,7 +1090,7 @@ subroutine AC_remesh_limitO2_pter(pos_adim, sca, bl_type, limit, buffer)
end if
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
bM=0.5*((y0-0.5)**2) - limit(1)
@@ -1142,10 +1142,10 @@ subroutine AC_remesh_limitO2_tag_CL_array(dir, pos_adim, sca, posP_ad, scaP, lim
j0_bis = floor(posP_ad)
!j0_bis = floor(posP/d_sc(2))
- y0 = (pos_adim - dble(j0))
- !y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
- y0_bis = (posP_ad - dble(j0_bis))
- !y0_bis = (posP - dble(j0_bis)*d_sc(2))/d_sc(2)
+ y0 = (pos_adim - real(j0, WP))
+ !y0 = (pos - real(j0, WP)*d_sc(2))/d_sc(2)
+ y0_bis = (posP_ad - real(j0_bis, WP))
+ !y0_bis = (posP - real(j0_bis, WP)*d_sc(2))/d_sc(2)
aM=0.5*((y0-0.5)**2) - limit(1) ! = (lambda 2 limited) alpha(y0_bis)
a0=1.-aM
@@ -1201,10 +1201,10 @@ subroutine AC_remesh_limitO2_tag_CL_pter(pos_adim, sca, posP_ad, scaP, limit, bu
jM=j0-1
jP=j0+1
- y0 = (pos_adim - dble(j0))
- !y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
- y0_bis = (posP_ad - dble(j0_bis))
- !y0_bis = (posP - dble(j0_bis)*d_sc(2))/d_sc(2)
+ y0 = (pos_adim - real(j0, WP))
+ !y0 = (pos - real(j0, WP)*d_sc(2))/d_sc(2)
+ y0_bis = (posP_ad - real(j0_bis, WP))
+ !y0_bis = (posP - real(j0_bis, WP)*d_sc(2))/d_sc(2)
aM=0.5*((y0-1.)**2) - limit(1) ! = (lambda 2 limited) alpha(y0_bis)
a0=1.-aM
@@ -1264,10 +1264,10 @@ subroutine AC_remesh_limitO2_tag_LC_array(dir, pos_adim, sca, posP_ad, scaP, lim
jP3=j0+3
! Distance to mesh point
- y0 = (pos_adim - dble(j0))
- !y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
- y0_bis = (posP_ad - dble(j0_bis))
- !y0_bis = (posP - dble(j0_bis)*d_sc(2))/d_sc(2)
+ y0 = (pos_adim - real(j0, WP))
+ !y0 = (pos - real(j0, WP)*d_sc(2))/d_sc(2)
+ y0_bis = (posP_ad - real(j0_bis, WP))
+ !y0_bis = (posP - real(j0_bis, WP)*d_sc(2))/d_sc(2)
! Interpolation weight
! Use limit(1) and limit(2) to remesh particle i (they are limitator at i-1/2, i+1/2)
@@ -1338,10 +1338,10 @@ subroutine AC_remesh_limitO2_tag_LC_pter(pos_adim, sca, posP_ad, scaP, limit, bu
jP3=j0+3
! Distance to mesh point
- y0 = (pos_adim - dble(j0))
- !y0 = (pos - dble(j0)*d_sc(2))/d_sc(2)
- y0_bis = (posP_ad - dble(j0_bis))
- !y0_bis = (posP - dble(j0_bis)*d_sc(2))/d_sc(2)
+ y0 = (pos_adim - real(j0, WP))
+ !y0 = (pos - real(j0, WP)*d_sc(2))/d_sc(2)
+ y0_bis = (posP_ad - real(j0_bis, WP))
+ !y0_bis = (posP - real(j0_bis, WP)*d_sc(2))/d_sc(2)
! Interpolation weight
! Use limit(1) and limit(2) to remesh particle i (they are limitator at i-1/2, i+1/2)
@@ -1392,7 +1392,7 @@ subroutine AC_remesh_O4_left_array(dir, pos_adim, sca, buffer)
!j0 = floor(pos/d_sc(2))
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
!bM2=(y0-2.)*(y0-1.)*y0*(y0+1.)/24.0
@@ -1428,7 +1428,7 @@ subroutine AC_remesh_O4_left_pter(pos_adim, sca, buffer)
use cart_topology
use advec_variables ! contains info about solver parameters and others.
-
+
!Input/Ouput
real(WP), intent(in) :: pos_adim, sca
type(real_pter), dimension(:), intent(inout) :: buffer
@@ -1439,9 +1439,9 @@ subroutine AC_remesh_O4_left_pter(pos_adim, sca, buffer)
! Mesh point used in remeshing formula
j0 = floor(pos_adim)
!j0 = floor(pos/d_sc(2))
-
+
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
!bM2=(y0-2.)*(y0-1.)*y0*(y0+1.)/24.0
@@ -1486,7 +1486,7 @@ subroutine AC_remesh_O4_center_array(dir, pos_adim, sca, buffer)
j0 = nint(pos_adim)
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
!bM2=(y0-2.)*(y0-1.)*y0*(y0+1.)/24.0
@@ -1532,9 +1532,9 @@ subroutine AC_remesh_O4_center_pter(pos_adim, sca, buffer)
real(WP) :: y0 ! adimensionned distance to mesh points
! Mesh point used in remeshing formula
j0 = nint(pos_adim)
-
+
! Distance to mesh points
- y0 = (pos_adim - dble(j0))
+ y0 = (pos_adim - real(j0, WP))
! Interpolation weights
!bM2=(y0-2.)*(y0-1.)*y0*(y0+1.)/24.0
@@ -1555,7 +1555,7 @@ subroutine AC_remesh_O4_center_pter(pos_adim, sca, buffer)
buffer(j0+1)%pter = buffer(j0+1)%pter + bP*sca
buffer(j0+2)%pter = buffer(j0+2)%pter + bP2*sca
-
+
end subroutine AC_remesh_O4_center_pter
@@ -1603,10 +1603,10 @@ subroutine AC_remesh_O4_tag_CL_array(dir, posM_ad, scaM, pos_adim, sca, posP_ad,
jP2= floor(posP2_ad)
! Distance to mesh point
- yM = (posM_ad - dble(jM))
- y0 = (pos_adim - dble(j0))
- yP = (posP_ad - dble(jP))
- yP2= (posP2_ad - dble(jP2))
+ yM = (posM_ad - real(jM, WP))
+ y0 = (pos_adim - real(j0, WP))
+ yP = (posP_ad - real(jP, WP))
+ yP2= (posP2_ad - real(jP2, WP))
! Interpolation weights
!aM3=(yM-2.)*(yM-1.)*yM*(yM+1.)/24.0
@@ -1622,8 +1622,8 @@ subroutine AC_remesh_O4_tag_CL_array(dir, posM_ad, scaM, pos_adim, sca, posP_ad,
bM2=y0*(2.+y0*(-1.+y0*(-2.+y0)))/24.0
!bM =(2.-y0)*(y0-1.)*y0*(y0+2.)/6.0
bM =y0*(-4.+y0*(4.+y0*(1.-y0)))/6.0
- !bP =((y0+1)-1.)*(y0+1)*((y0+1)+1.)*((y0+1)+2.)/24.0
- bP =y0*(6.+y0*(11+y0*(6+y0)))/24.0
+ !bP =((y0+1)-1.)*(y0+1)*((y0+1)+1.)*((y0+1)+2.)/24.0
+ bP =y0*(6.+y0*(11+y0*(6+y0)))/24.0
!b0 =((y0-2.)*(y0-1.)*(y0+1.)*(y0+2.)/4.0) + ((2.-y0)*y0*(y0+1.)*(y0+2.)/6.0) &
! & + ((y0-1.)*y0*(y0+1.)*(y0+2.)/24.0) - bP
b0 = 1. - (bM2+bM+bP)
@@ -1647,12 +1647,12 @@ subroutine AC_remesh_O4_tag_CL_array(dir, posM_ad, scaM, pos_adim, sca, posP_ad,
!e0 =((yP2-2.)*(yP2-1.)*yP2*(yP2+1.)/24.0) + ((2.-yP2)*(yP2-1.)*yP2*(yP2+2.)/6.0)
e0 = 1. - (eP+eP2+eP3)
- ! -- remeshing --
+ ! -- remeshing --
! j0-3
j1 = modulo(j0-4,mesh_sc%N(dir))+1
buffer(j1) = buffer(j1) + aM3*scaM
! j0-2
- j1 = modulo(j0-3,mesh_sc%N(dir))+1
+ j1 = modulo(j0-3,mesh_sc%N(dir))+1
buffer(j1) = buffer(j1) + aM2*scaM + bM2*sca
! j0-1
j1 = modulo(j0-2,mesh_sc%N(dir))+1
@@ -1715,10 +1715,10 @@ subroutine AC_remesh_O4_tag_CL_pter(posM_ad, scaM, pos_adim, sca, posP_ad, scaP,
jP2= floor(posP2_ad)
! Distance to mesh point
- yM = (posM_ad - dble(jM))
- y0 = (pos_adim - dble(j0))
- yP = (posP_ad - dble(jP))
- yP2= (posP2_ad - dble(jP2))
+ yM = (posM_ad - real(jM, WP))
+ y0 = (pos_adim - real(j0, WP))
+ yP = (posP_ad - real(jP, WP))
+ yP2= (posP2_ad - real(jP2, WP))
! Interpolation weights
!aM3=(yM-2.)*(yM-1.)*yM*(yM+1.)/24.0
@@ -1734,7 +1734,7 @@ subroutine AC_remesh_O4_tag_CL_pter(posM_ad, scaM, pos_adim, sca, posP_ad, scaP,
bM2=y0*(2._WP+y0*(-1._WP+y0*(-2._WP+y0)))/24._WP
!bM =(2.-y0)*(y0-1.)*y0*(y0+2.)/6.0
bM =y0*(-4._WP+y0*(4._WP+y0*(1._WP-y0)))/6._WP
- !bP =((y0+1)-1.)*(y0+1)*((y0+1)+1.)*((y0+1)+2.)/24.0
+ !bP =((y0+1)-1.)*(y0+1)*((y0+1)+1.)*((y0+1)+2.)/24.0
bP =y0*(6._WP+y0*(11._WP+y0*(6._WP+y0)))/24._WP
!b0 =((y0-2.)*(y0-1.)*(y0+1.)*(y0+2.)/4.0) + ((2.-y0)*y0*(y0+1.)*(y0+2.)/6.0) &
! & + ((y0-1.)*y0*(y0+1.)*(y0+2.)/24.0) - bP
@@ -1760,7 +1760,7 @@ subroutine AC_remesh_O4_tag_CL_pter(posM_ad, scaM, pos_adim, sca, posP_ad, scaP,
e0 = 1._WP - (eP+eP2+eP3)
! remeshing
- buffer(j0-3)%pter = buffer(j0-3)%pter +aM3*scaM
+ buffer(j0-3)%pter = buffer(j0-3)%pter +aM3*scaM
buffer(j0-2)%pter = buffer(j0-2)%pter +aM2*scaM +bM2*sca
buffer(j0-1)%pter = buffer(j0-1)%pter + aM*scaM + bM*sca + cM*scaP
buffer(j0 )%pter = buffer(j0 )%pter + a0*scaM + b0*sca + c0*scaP + e0*scaP2
@@ -1818,10 +1818,10 @@ subroutine AC_remesh_O4_tag_LC_array(dir, posM_ad, scaM, pos_adim, sca, posP_ad,
jP2= nint(posP2_ad)
! Distance to mesh point
- yM = (posM_ad - dble(jM))
- y0 = (pos_adim - dble(j0))
- yP = (posP_ad - dble(jP))
- yP2= (posP2_ad - dble(jP2))
+ yM = (posM_ad - real(jM, WP))
+ y0 = (pos_adim - real(j0, WP))
+ yP = (posP_ad - real(jP, WP))
+ yP2= (posP2_ad - real(jP2, WP))
! Interpolation weights
!aM3=(yM-2.)*(yM-1.)*yM*(yM+1.)/24.0
@@ -1830,7 +1830,7 @@ subroutine AC_remesh_O4_tag_LC_array(dir, posM_ad, scaM, pos_adim, sca, posP_ad,
aM2 =yM*(-4.+yM*(4.+yM*(1.-yM)))/6.0
!aM =(yM-2.)*(yM-1.)*(yM+1.)*(yM+2.)/4.0
aM =(4.+(yM**2)*(-5.+yM**2))/4.0
- !a0 =((2.-yM)*yM*(yM+1.)*(yM+2.)/6.0)
+ !a0 =((2.-yM)*yM*(yM+1.)*(yM+2.)/6.0)
a0 =yM*(4+yM*(4-yM*(1.+yM)))/6.0
!aP2=(((yM-1.)-1.)*(yM-1.)*((yM-1.)+1.)*((yM-1.)+2.)/24.0)
!aP2=yM*(yM-2.)*(yM-1.)*(yM+1.)/24.0
@@ -1962,10 +1962,10 @@ subroutine AC_remesh_O4_tag_LC_pter(posM_ad, scaM, pos_adim, sca, posP_ad, scaP,
jP2= nint(posP2_ad)
! Distance to mesh point
- yM = (posM_ad - dble(jM))
- y0 = (pos_adim - dble(j0))
- yP = (posP_ad - dble(jP))
- yP2= (posP2_ad - dble(jP2))
+ yM = (posM_ad - real(jM, WP))
+ y0 = (pos_adim - real(j0, WP))
+ yP = (posP_ad - real(jP, WP))
+ yP2= (posP2_ad - real(jP2, WP))
! Interpolation weights
!aM3=(yM-2.)*(yM-1.)*yM*(yM+1.)/24.0
@@ -1974,7 +1974,7 @@ subroutine AC_remesh_O4_tag_LC_pter(posM_ad, scaM, pos_adim, sca, posP_ad, scaP,
aM2 =yM*(-4._WP+yM*(4._WP+yM*(1._WP-yM)))/6._WP
!aM =(yM-2.)*(yM-1.)*(yM+1.)*(yM+2.)/4.0
aM =(4._WP+(yM**2)*(-5._WP+yM**2))/4._WP
- !a0 =((2.-yM)*yM*(yM+1.)*(yM+2.)/6.0)
+ !a0 =((2.-yM)*yM*(yM+1.)*(yM+2.)/6.0)
a0 =yM*(4._WP+yM*(4._WP-yM*(1._WP+yM)))/6._WP
!aP2=(((yM-1.)-1.)*(yM-1.)*((yM-1.)+1.)*((yM-1.)+2.)/24.0)
!aP2=yM*(yM-2.)*(yM-1.)*(yM+1.)/24.0
diff --git a/HySoP/src/scalesInterface/particles/interpolation_velo.f90 b/HySoP/src/scalesInterface/particles/interpolation_velo.f90
index 51290ff88..3acada44e 100644
--- a/HySoP/src/scalesInterface/particles/interpolation_velo.f90
+++ b/HySoP/src/scalesInterface/particles/interpolation_velo.f90
@@ -453,7 +453,7 @@ subroutine Inter_LastDir_com(dir, V_coarse, dx_c, V_fine, dx_f)
! ==== Communication ====
if(stencil_g>0) then
allocate(V_beg(size(V_coarse,1),size(V_coarse,2),stencil_g))
- com_nb(1) = ceiling(real(stencil_g)/N_coarse) ! number of required communication to get ghost
+ com_nb(1) = ceiling(real(stencil_g, WP)/N_coarse) ! number of required communication to get ghost
allocate(beg_request(com_nb(1)))
com_pos = stencil_g+1 ! i = 1 + missing (or remainding) ghost lines
com_size(2) = com_size(1)*N_coarse
@@ -479,7 +479,7 @@ subroutine Inter_LastDir_com(dir, V_coarse, dx_c, V_fine, dx_f)
if(stencil_d>0) then
allocate(V_end(size(V_coarse,1),size(V_coarse,2),stencil_d))
- com_nb(2) = ceiling(real(stencil_d)/N_coarse) ! number of required communication to get ghost
+ com_nb(2) = ceiling(real(stencil_d, WP)/N_coarse) ! number of required communication to get ghost
allocate(end_request(com_nb(2)))
com_pos = 1 ! Reception from next processus is done in position 1
com_size(2) = com_size(1)*N_coarse
@@ -625,7 +625,7 @@ subroutine Inter_LastDir_Permut_com(dir, V_coarse, dx_c, V_fine, dx_f)
! ==== Communication ====
if(stencil_g>0) then
allocate(V_beg(size(V_coarse,1),size(V_coarse,2),stencil_g))
- com_nb(1) = ceiling(real(stencil_g)/N_coarse) ! number of required communication to get ghost
+ com_nb(1) = ceiling(real(stencil_g, WP)/N_coarse) ! number of required communication to get ghost
allocate(beg_request(com_nb(1)))
com_pos = stencil_g+1 ! i = 1 + missing (or remainding) ghost lines
com_size(2) = com_size(1)*N_coarse
@@ -651,7 +651,7 @@ subroutine Inter_LastDir_Permut_com(dir, V_coarse, dx_c, V_fine, dx_f)
if(stencil_d>0) then
allocate(V_end(size(V_coarse,1),size(V_coarse,2),stencil_d))
- com_nb(2) = ceiling(real(stencil_d)/N_coarse) ! number of required communication to get ghost
+ com_nb(2) = ceiling(real(stencil_d, WP)/N_coarse) ! number of required communication to get ghost
allocate(end_request(com_nb(2)))
com_pos = 1 ! Reception from next processus is done in position 1
com_size(2) = com_size(1)*N_coarse
diff --git a/HySoP/src/scalesInterface/precision_tools.f90 b/HySoP/src/scalesInterface/precision_tools.f90
index 2818f9ad9..9e6213968 100644
--- a/HySoP/src/scalesInterface/precision_tools.f90
+++ b/HySoP/src/scalesInterface/precision_tools.f90
@@ -34,7 +34,7 @@ MODULE precision_tools
INTEGER, PUBLIC :: MPI_REAL_WP
!> the MPI type for COMPLEX exchanges in simple or double precision
INTEGER, PUBLIC :: MPI_COMPLEX_WP
- !> the string size
+ !> the string size
INTEGER, PARAMETER :: str_short = 8
INTEGER, PARAMETER :: str_medium = 64
INTEGER, PARAMETER :: str_long = 4096
--
GitLab