From 9736fd28dd4c533928c9597d9860d8fc0609b60b Mon Sep 17 00:00:00 2001
From: jtra <jtra@7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9>
Date: Tue, 19 Apr 2011 10:46:47 +0000
Subject: [PATCH] Included new fftw wrapper routines Included poisson solver
that depends on the new fftw wrappers
git-svn-id: https://ppm.inf.ethz.ch/svn/ppmnumerics/branches/ngtopo/libppmnumerics@840 7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9
---
src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f90 | 117 ++++
src/fft/ppm_fft_exec_3d_vec_c2c_z.f90 | 122 ++++
src/fft/ppm_fft_exec_3d_vec_fr2c_xy.f90 | 103 ++++
src/fft/ppm_fft_normalize_c.f90 | 119 ++++
src/fft/ppm_fft_normalize_r.f90 | 120 ++++
src/fft/ppm_fft_plan_3d_vec_bc2c_z.f90 | 145 +++++
src/fft/ppm_fft_plan_3d_vec_bc2r_xy.f90 | 149 +++++
src/fft/ppm_fft_plan_3d_vec_fc2c_z.f90 | 145 +++++
src/fft/ppm_fft_plan_3d_vec_fr2c_xy.f90 | 148 +++++
src/poisson/ppm_poisson_fd.f90 | 144 +++++
src/poisson/ppm_poisson_init_predef.f90 | 414 ++++++++++++++
src/poisson/ppm_poisson_solve.f90 | 712 ++++++++++++++++++++++++
src/ppm_define.h | 2 +
src/ppm_module_fft.f90 | 154 +++++
src/ppm_module_poisson.f90 | 195 +++++++
15 files changed, 2789 insertions(+)
create mode 100644 src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f90
create mode 100644 src/fft/ppm_fft_exec_3d_vec_c2c_z.f90
create mode 100644 src/fft/ppm_fft_exec_3d_vec_fr2c_xy.f90
create mode 100644 src/fft/ppm_fft_normalize_c.f90
create mode 100644 src/fft/ppm_fft_normalize_r.f90
create mode 100644 src/fft/ppm_fft_plan_3d_vec_bc2c_z.f90
create mode 100644 src/fft/ppm_fft_plan_3d_vec_bc2r_xy.f90
create mode 100644 src/fft/ppm_fft_plan_3d_vec_fc2c_z.f90
create mode 100644 src/fft/ppm_fft_plan_3d_vec_fr2c_xy.f90
create mode 100644 src/poisson/ppm_poisson_fd.f90
create mode 100644 src/poisson/ppm_poisson_init_predef.f90
create mode 100644 src/poisson/ppm_poisson_solve.f90
create mode 100644 src/ppm_module_fft.f90
create mode 100644 src/ppm_module_poisson.f90
diff --git a/src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f90 b/src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f90
new file mode 100644
index 0000000..427f434
--- /dev/null
+++ b/src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f90
@@ -0,0 +1,117 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_fft_exec_3d_vec_bc2r_xy
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! FFTW execute wrapper for 3d arrays, 2d complex to real
+ ! (backward) FFT in the xy directions
+ ! The routine does not work with fields that include ghost layers
+ !-------------------------------------------------------------------------
+#if __KIND == __SINGLE
+#define __ROUTINE ppm_fft_exec_3d_vec_bc2r_xy_s
+#define __PREC ppm_kind_single
+#elif __KIND == __DOUBLE
+#define __ROUTINE ppm_fft_exec_3d_vec_bc2r_xy_d
+#define __PREC ppm_kind_double
+#endif
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmplan,infield,outfield,info)
+ !!! FFTW execute wrapper for 3d arrays, 2d complex to real
+ !!! (backward) FFT in the xy directions
+ !!! Before calling this routine a ppm_fft_plan_ routine must be called
+ !!! The routine does not work with fields that include ghost layers
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_typedef
+ USE ppm_module_topo_get
+ USE ppm_module_write
+ USE ppm_module_data,ONLY:ppm_rank,ppm_kind_single,ppm_kind_double
+
+ IMPLICIT NONE
+
+ INCLUDE 'fftw3.f'
+
+ ! if debug check if dimensions are 2a 3b 5c 7d 11e 13f
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ !!!topology identifier of target
+ INTEGER,INTENT(IN) :: topoid
+ !!!id of the mesh
+ INTEGER,INTENT(IN) :: meshid
+ !!!ppm fft plan type
+ TYPE(ppm_fft_plan),INTENT(INOUT) :: ppmplan
+ !!!input field to fourier transform
+ !COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: infield
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: infield
+ !!!output field for the result of the fourier transform
+ !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: outfield
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: outfield
+ !!!Returns status, 0 upon success
+ INTEGER,INTENT(OUT) :: info
+ !in time perhaps an argument for alternate directions
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ REAL(__PREC) :: t0
+ INTEGER :: i,k
+ INTEGER :: isub,isubl
+ INTEGER :: nsubs
+ INTEGER,DIMENSION(:),POINTER :: isublist
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_fft_exec',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_exec','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ nsubs = topology%nsublist
+ ALLOCATE(isublist(nsubs))
+ DO isub=1,nsubs
+ isublist(isub) = topology%isublist(isub)
+ ENDDO
+ mesh = topology%mesh(meshid)
+
+ !-------------------------------------------------------------------------
+ ! Execute plan
+ !-------------------------------------------------------------------------
+ DO isub=1,nsubs
+ isubl=isublist(isub)
+ DO k=1,mesh%nnodes(3,isubl) !@ add '-1' to exclude n+1 slabs
+ CALL dfftw_execute_dft_c2r(ppmplan%plan(isub),&
+ & infield(1,1,1,k,isub),outfield(1,1,1,k,isub))
+ !-------------------------------------------------------------------
+ ! Copy periodic layer back - only for periodic 'N' vertex points
+ !-------------------------------------------------------------------
+ DO i=1,mesh%nnodes(1,isubl)
+ outfield(1,i,mesh%nnodes(2,isubl),k,isub) = outfield(1,i,1,k,isub)
+ outfield(2,i,mesh%nnodes(2,isubl),k,isub) = outfield(2,i,1,k,isub)
+ outfield(3,i,mesh%nnodes(2,isubl),k,isub) = outfield(3,i,1,k,isub)
+ END DO
+ DO i=1,mesh%nnodes(2,isubl)
+ outfield(1,mesh%nnodes(1,isubl),i,k,isub) = outfield(1,1,i,k,isub)
+ outfield(2,mesh%nnodes(1,isubl),i,k,isub) = outfield(2,1,i,k,isub)
+ outfield(3,mesh%nnodes(1,isubl),i,k,isub) = outfield(3,1,i,k,isub)
+ END DO
+ END DO
+ END DO
+
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_fft_exec',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
diff --git a/src/fft/ppm_fft_exec_3d_vec_c2c_z.f90 b/src/fft/ppm_fft_exec_3d_vec_c2c_z.f90
new file mode 100644
index 0000000..f307aab
--- /dev/null
+++ b/src/fft/ppm_fft_exec_3d_vec_c2c_z.f90
@@ -0,0 +1,122 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_fft_exec_3d_vec_c2c_z
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! FFTW execute wrapper for 3d arrays, 1d complex to complex
+ ! (forward and backward) FFT in the z direction
+ ! The routine does not work with fields that include ghost layers
+ !-------------------------------------------------------------------------
+#if __KIND == __SINGLE
+#define __ROUTINE ppm_fft_exec_3d_vec_c2c_z_s
+#define __PREC ppm_kind_single
+#elif __KIND == __DOUBLE
+#define __ROUTINE ppm_fft_exec_3d_vec_c2c_z_d
+#define __PREC ppm_kind_double
+#endif
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmplan,infield,outfield,info)
+ !!! FFTW execute wrapper for 3d arrays, 1d complex to complex
+ !!! (forward and backward) FFT in the z direction
+ !!! Before calling this routine a ppm_fft_plan_ routine must be called
+ !!! The routine does not work with fields that include ghost layers
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_typedef
+ USE ppm_module_topo_get
+ USE ppm_module_write
+ USE ppm_module_data,ONLY:ppm_rank,ppm_kind_single,ppm_kind_double
+
+ IMPLICIT NONE
+
+ INCLUDE 'fftw3.f'
+
+ ! if debug check if dimensions are 2a 3b 5c 7d 11e 13f
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ !!!topology identifier of target
+ INTEGER,INTENT(IN) :: topoid
+ !!!id of the mesh
+ INTEGER,INTENT(IN) :: meshid
+ !!!ppm fft plan type
+ TYPE(ppm_fft_plan),INTENT(INOUT) :: ppmplan
+ !!!input field to fourier transform
+ !COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: infield
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: infield
+ !!!output field for the result of the fourier transform
+ !COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: outfield
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: outfield
+ !!!Returns status, 0 upon success
+ INTEGER,INTENT(OUT) :: info
+ !in time perhaps an argument for alternate directions
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ REAL(__PREC) :: t0
+ INTEGER :: i,j
+ INTEGER :: isub,isubl
+ INTEGER :: nsubs
+ INTEGER,DIMENSION(:),POINTER :: isublist
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_fft_exec',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_exec','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ nsubs = topology%nsublist
+ ALLOCATE(isublist(nsubs))
+ DO isub=1,nsubs
+ isublist(isub) = topology%isublist(isub)
+ ENDDO
+ mesh = topology%mesh(meshid)
+
+ !outfield = 0.0_ppm_kind_double !@
+ !-------------------------------------------------------------------------
+ ! Execute plan
+ !-------------------------------------------------------------------------
+ DO isub=1,nsubs
+ isubl=isublist(isub)
+ DO j=1,mesh%nnodes(2,isubl) !@ add '-1' to exclude n+1 points
+ DO i=1,mesh%nnodes(1,isubl) !@ add '-1' to exclude n+1 points
+ CALL dfftw_execute_dft(ppmplan%plan(isub),&
+ & infield(1,i,j,1,isub),outfield(1,i,j,1,isub))
+ END DO
+ END DO
+ END DO
+
+ !-------------------------------------------------------------------------
+ ! Copy periodic - this is only for the periodic 'N+1' vertex points
+ !-------------------------------------------------------------------------
+ IF (ppmplan%sign .EQ. FFTW_BACKWARD) THEN
+ DO isub=1,nsubs
+ isubl=isublist(isub)
+ DO j=1,mesh%nnodes(2,isubl) !@ add '-1' to exclude n+1 points
+ DO i=1,mesh%nnodes(1,isubl) !@ add '-1' to exclude n+1 points
+ outfield(1,i,j,mesh%nnodes(3,isubl),isub) = outfield(1,i,j,1,isub)
+ outfield(2,i,j,mesh%nnodes(3,isubl),isub) = outfield(2,i,j,1,isub)
+ outfield(3,i,j,mesh%nnodes(3,isubl),isub) = outfield(3,i,j,1,isub)
+ END DO
+ END DO
+ END DO
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_fft_exec',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
diff --git a/src/fft/ppm_fft_exec_3d_vec_fr2c_xy.f90 b/src/fft/ppm_fft_exec_3d_vec_fr2c_xy.f90
new file mode 100644
index 0000000..82ea0e3
--- /dev/null
+++ b/src/fft/ppm_fft_exec_3d_vec_fr2c_xy.f90
@@ -0,0 +1,103 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_fft_exec_3d_vec_fr2c_xy
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! FFTW execute wrapper for 3d arrays, 2d real to complex
+ ! (forward) FFT in the xy directions
+ ! The routine does not work with fields that include ghost layers
+ !-------------------------------------------------------------------------
+#if __KIND == __SINGLE
+#define __ROUTINE ppm_fft_exec_3d_vec_fr2c_xy_s
+#define __PREC ppm_kind_single
+#elif __KIND == __DOUBLE
+#define __ROUTINE ppm_fft_exec_3d_vec_fr2c_xy_d
+#define __PREC ppm_kind_double
+#endif
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmplan,infield,outfield,info)
+ !!! FFTW execute wrapper for 3d arrays, 2d real to complex
+ !!! (forward) FFT in the xy directions
+ !!! Before calling this routine a ppm_fft_plan_ routine must be called
+ !!! The routine does not work with fields that include ghost layers
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_typedef
+ USE ppm_module_topo_get
+ USE ppm_module_write
+ USE ppm_module_data,ONLY:ppm_rank,ppm_kind_single,ppm_kind_double
+
+ IMPLICIT NONE
+
+ INCLUDE 'fftw3.f'
+
+ ! if debug check if dimensions are 2a 3b 5c 7d 11e 13f
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ !!!topology identifier of target
+ INTEGER,INTENT(IN) :: topoid
+ !!!id of the mesh
+ INTEGER,INTENT(IN) :: meshid
+ !!!ppm fft plan type
+ TYPE(ppm_fft_plan),INTENT(INOUT) :: ppmplan
+ !!!input field to fourier transform
+ !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: infield
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: infield
+ !!!output field for the result of the fourier transform
+ !COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: outfield
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: outfield
+ !!!Returns status, 0 upon success
+ INTEGER,INTENT(OUT) :: info
+ !in time perhaps an argument for alternate directions
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ REAL(__PREC) :: t0
+ INTEGER :: k
+ INTEGER :: isub,isubl
+ INTEGER :: nsubs
+ INTEGER,DIMENSION(:),POINTER :: isublist
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_fft_exec',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_exec','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ nsubs = topology%nsublist
+ ALLOCATE(isublist(nsubs))
+ DO isub=1,nsubs
+ isublist(isub) = topology%isublist(isub)
+ ENDDO
+ mesh = topology%mesh(meshid)
+
+ !-------------------------------------------------------------------------
+ ! Execute plan
+ !-------------------------------------------------------------------------
+ DO isub=1,nsubs
+ isubl=isublist(isub)
+ DO k=1,mesh%nnodes(3,isubl) !@ add '-1' to exclude n+1 slabs
+ CALL dfftw_execute_dft_r2c(ppmplan%plan(isub),&
+ & infield(1,1,1,k,isub),outfield(1,1,1,k,isub))
+ END DO
+ END DO
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_fft_exec',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
diff --git a/src/fft/ppm_fft_normalize_c.f90 b/src/fft/ppm_fft_normalize_c.f90
new file mode 100644
index 0000000..a9d18bb
--- /dev/null
+++ b/src/fft/ppm_fft_normalize_c.f90
@@ -0,0 +1,119 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_fft_normalize_c
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! The routine does not work with fields that include ghost layers.
+ ! Fields are noralized in the 1/N - fashion that scales fields to the
+ ! proper level after being FFTed and iFFTed
+ ! The routines are not meant for production code as the normalization can
+ ! be done for all dimensions in one pass - which may and should be done
+ ! in a routine (looping through the data anyway) following the FFTs.
+ !-------------------------------------------------------------------------
+#if __KIND == __SINGLE
+#define __ROUTINE ppm_fft_normalize_cs
+#define __PREC ppm_kind_single
+#elif __KIND == __DOUBLE
+#define __ROUTINE ppm_fft_normalize_cd
+#define __PREC ppm_kind_double
+#endif
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmplan,infield,info)
+ !!! The routine does not work with fields that include ghost layers.
+ !!! Fields are noralized in the 1/N - fashion that scales fields to the
+ !!! proper level after being FFTed and iFFTed
+ !!! The routines are not meant for production code as the normalization
+ !!! can be done for all dimensions in one pass - which may and should be
+ !!! done in a routine (looping through the data anyway) following the FFTs
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_typedef
+ USE ppm_module_topo_get
+ USE ppm_module_write
+ USE ppm_module_data,ONLY:ppm_rank,ppm_kind_single,ppm_kind_double
+
+ IMPLICIT NONE
+
+ INCLUDE 'fftw3.f'
+
+ ! if debug check if dimensions are 2a 3b 5c 7d 11e 13f
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ !!!topology identifier of target
+ INTEGER,INTENT(IN) :: topoid
+ !!!id of the mesh
+ INTEGER,INTENT(IN) :: meshid
+ !!!ppm fft plan type
+ TYPE(ppm_fft_plan),INTENT(INOUT) :: ppmplan
+ !!!input field to normalize
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: infield
+ !!!Returns status, 0 upon success
+ INTEGER,INTENT(OUT) :: info
+ !in time perhaps an argument for alternate directions
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ !timer
+ INTEGER,PARAMETER :: MK = ppm_kind_double
+ REAL(__PREC) :: t0
+ !normalization factor
+ REAL(MK) :: fac
+ INTEGER :: i,j,k
+ INTEGER :: isub,isubl
+ INTEGER :: nsubs
+ INTEGER,DIMENSION(:),POINTER :: isublist
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_fft_normalize',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_plan','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ nsubs = topology%nsublist
+ ALLOCATE(isublist(nsubs))
+ DO isub=1,nsubs
+ isublist(isub) = topology%isublist(isub)
+ ENDDO
+ mesh = topology%mesh(meshid)
+
+ DO isub=1,nsubs
+ isubl=isublist(isub)
+ !determine normalization factor
+ IF (ppmplan%rank .EQ. 1) THEN
+ fac = 1.0_MK/REAL((mesh%nm(3)-1),MK)
+ ELSE IF (ppmplan%rank .EQ. 2) THEN
+ fac = 1.0_MK/REAL((mesh%nm(1)-1)*(mesh%nm(2)-1),MK)
+ ENDIF
+
+ DO k=1,mesh%nnodes(3,isubl)
+ DO j=1,mesh%nnodes(2,isubl)
+ DO i=1,mesh%nnodes(1,isubl)
+ infield(1,i,j,k,isub) = fac*infield(1,i,j,k,isub)
+ infield(2,i,j,k,isub) = fac*infield(2,i,j,k,isub)
+ infield(3,i,j,k,isub) = fac*infield(3,i,j,k,isub)
+ END DO
+ END DO
+ END DO
+ END DO
+
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_fft_normalize',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
+
diff --git a/src/fft/ppm_fft_normalize_r.f90 b/src/fft/ppm_fft_normalize_r.f90
new file mode 100644
index 0000000..0591d33
--- /dev/null
+++ b/src/fft/ppm_fft_normalize_r.f90
@@ -0,0 +1,120 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_fft_normalize_r
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! The routine does not work with fields that include ghost layers.
+ ! Fields are noralized in the 1/N - fashion that scales fields to the
+ ! proper level after being FFTed and iFFTed
+ ! The routines are not meant for production code as the normalization can
+ ! be done for all dimensions in one pass - which may and should be done
+ ! in a routine (looping through the data anyway) following the FFTs.
+ !-------------------------------------------------------------------------
+#if __KIND == __SINGLE
+#define __ROUTINE ppm_fft_normalize_rs
+#define __PREC ppm_kind_single
+#elif __KIND == __DOUBLE
+#define __ROUTINE ppm_fft_normalize_rd
+#define __PREC ppm_kind_double
+#endif
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmplan,infield,info)
+ !!! The routine does not work with fields that include ghost layers.
+ !!! Fields are noralized in the 1/N - fashion that scales fields to the
+ !!! proper level after being FFTed and iFFTed
+ !!! The routines are not meant for production code as the normalization
+ !!! can be done for all dimensions in one pass - which may and should be
+ !!! done in a routine (looping through the data anyway) following the FFTs
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_typedef
+ USE ppm_module_topo_get
+ USE ppm_module_write
+ USE ppm_module_data,ONLY:ppm_rank,ppm_kind_single,ppm_kind_double
+
+ IMPLICIT NONE
+
+ INCLUDE 'fftw3.f'
+
+ ! if debug check if dimensions are 2a 3b 5c 7d 11e 13f
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ !!!topology identifier of target
+ INTEGER,INTENT(IN) :: topoid
+ !!!id of the mesh
+ INTEGER,INTENT(IN) :: meshid
+ !!!ppm fft plan type
+ TYPE(ppm_fft_plan),INTENT(INOUT) :: ppmplan
+ !!!input field to normalize
+ !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: infield
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: infield
+ !!!Returns status, 0 upon success
+ INTEGER,INTENT(OUT) :: info
+ !in time perhaps an argument for alternate directions
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ !timer
+ INTEGER,PARAMETER :: MK = ppm_kind_double
+ REAL(__PREC) :: t0
+ !normalization factor
+ REAL(MK) :: fac
+ INTEGER :: i,j,k
+ INTEGER :: isub,isubl
+ INTEGER :: nsubs
+ INTEGER,DIMENSION(:),POINTER :: isublist
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_fft_normalize',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_plan','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ nsubs = topology%nsublist
+ ALLOCATE(isublist(nsubs))
+ DO isub=1,nsubs
+ isublist(isub) = topology%isublist(isub)
+ ENDDO
+ mesh = topology%mesh(meshid)
+
+ DO isub=1,nsubs
+ isubl=isublist(isub)
+ !determine normalization factor
+ IF (ppmplan%rank .EQ. 1) THEN
+ fac = 1.0_MK/REAL((mesh%nm(3)-1),MK)
+ ELSE IF (ppmplan%rank .EQ. 2) THEN
+ fac = 1.0_MK/REAL((mesh%nm(1)-1)*(mesh%nm(2)-1),MK)
+ ENDIF
+
+ DO k=1,mesh%nnodes(3,isubl)
+ DO j=1,mesh%nnodes(2,isubl)
+ DO i=1,mesh%nnodes(1,isubl)
+ infield(1,i,j,k,isub) = fac*infield(1,i,j,k,isub)
+ infield(2,i,j,k,isub) = fac*infield(2,i,j,k,isub)
+ infield(3,i,j,k,isub) = fac*infield(3,i,j,k,isub)
+ END DO
+ END DO
+ END DO
+ END DO
+
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_fft_normalize',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
+
diff --git a/src/fft/ppm_fft_plan_3d_vec_bc2c_z.f90 b/src/fft/ppm_fft_plan_3d_vec_bc2c_z.f90
new file mode 100644
index 0000000..c08e3a8
--- /dev/null
+++ b/src/fft/ppm_fft_plan_3d_vec_bc2c_z.f90
@@ -0,0 +1,145 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_fft_plan_3d_vec_bc2c_z
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! FFTW plan wrapper for 3d arrays, 1d complex to complex
+ ! (backward) FFT in the z direction
+ ! The routine does not work with fields that include ghost layers
+ !-------------------------------------------------------------------------
+#if __KIND == __SINGLE
+#define __ROUTINE ppm_fft_plan_3d_vec_bc2c_z_s
+#define __PREC ppm_kind_single
+#elif __KIND == __DOUBLE
+#define __ROUTINE ppm_fft_plan_3d_vec_bc2c_z_d
+#define __PREC ppm_kind_double
+#endif
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmplan,infield,outfield,info)
+ !!! FFTW plan wrapper for 3d arrays, 1d complex to complex
+ !!! (backward) FFT in the z direction
+ !!! The routine does not work with fields that include ghost layers
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_typedef
+ USE ppm_module_topo_get
+ USE ppm_module_write
+ USE ppm_module_data,ONLY:ppm_rank,ppm_kind_single,ppm_kind_double
+
+ IMPLICIT NONE
+
+ INCLUDE 'fftw3.f'
+
+ ! if debug check if dimensions are 2a 3b 5c 7d 11e 13f
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ !!!topology identifier of target
+ INTEGER,INTENT(IN) :: topoid
+ !!!id of the mesh
+ INTEGER,INTENT(IN) :: meshid
+ !!!ppm fft plan type
+ TYPE(ppm_fft_plan),INTENT(INOUT) :: ppmplan
+ !!!input field to fourier transform
+ !COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: infield
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: infield
+ !!!output field for the result of the fourier transform
+ !COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: outfield
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: outfield
+ !!!Returns status, 0 upon success
+ INTEGER,INTENT(OUT) :: info
+ !in time perhaps an argument for alternate directions
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ REAL(__PREC) :: t0
+ INTEGER :: isub,isubl
+ INTEGER :: nsubs
+ INTEGER,DIMENSION(:),POINTER :: isublist
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_fft_plan',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_plan','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ nsubs = topology%nsublist
+ ALLOCATE(isublist(nsubs))
+ DO isub=1,nsubs
+ isublist(isub) = topology%isublist(isub)
+ ENDDO
+ mesh = topology%mesh(meshid)
+
+ !-------------------------------------------------------------------------
+ ! Setup parameters for this particular routine
+ !-------------------------------------------------------------------------
+ !the dimension of the FFT (1D/2D/3D)
+ ppmplan%rank=1
+ !the number of points along each direction of the piece to be transformed
+ ALLOCATE(ppmplan%nx(ppmplan%rank,nsubs))
+ !the direction of the transform
+ ppmplan%sign=FFTW_BACKWARD
+ !the method to setup the optimal plan
+ ppmplan%flag=FFTW_MEASURE
+ !the number of components to transform - 3 component vector
+ ppmplan%howmany=3
+ !the size of the input array - full size (assuming LBOUND=1 thus UBOUND)
+ ALLOCATE(ppmplan%inembed(ppmplan%rank))
+ ppmplan%inembed(1) = UBOUND(infield,4)
+ !the size of the output array - full size (assuming LBOUND=1 thus UBOUND)
+ ALLOCATE(ppmplan%onembed(ppmplan%rank))
+ ppmplan%onembed(1) = UBOUND(outfield,4)
+ !istride tells how the same componenet data points are spaced in memory
+ !e.g. z values recur every x-dim*y-dim*component values
+ ppmplan%istride = UBOUND(infield,2) *UBOUND(infield,3)*3
+ ppmplan%ostride = UBOUND(outfield,2)*UBOUND(outfield,3)*3
+ !idist tells how multiple arrays are spaced in memory. I.e. a memory
+ !offset. e.g. vector components (idist=1) or scalar 2D arrays in
+ !3D array(idist=NxNy)
+ ppmplan%idist = 1
+ ppmplan%odist = 1
+
+ !-------------------------------------------------------------------------
+ ! Allocate plan array
+ !-------------------------------------------------------------------------
+ IF(ASSOCIATED(ppmplan%plan)) THEN
+ DEALLOCATE(ppmplan%plan,stat=info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_plan','Failed to deallocate plan-array.',isub)
+ GOTO 9999
+ ENDIF
+ END IF
+ ALLOCATE(ppmplan%plan(nsubs))
+
+ DO isub=1,nsubs
+ isubl=isublist(isub)
+ !@ maybe the -1 needs to be removed when doing cell data
+ !we subtract the -1 to avoid the periodic vertex point
+ ppmplan%nx(1,isub) = mesh%nnodes(3,isubl)-1
+
+ CALL dfftw_plan_many_dft(ppmplan%plan(isub),ppmplan%rank,&
+ & ppmplan%nx(:,isub),ppmplan%howmany,infield(1,1,1,1,isub),&
+ & ppmplan%inembed(1),ppmplan%istride,ppmplan%idist,&
+ & outfield(1,1,1,1,isub),ppmplan%onembed(1),ppmplan%ostride,&
+ & ppmplan%odist,ppmplan%sign,ppmplan%flag)
+ END DO
+
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_fft_plan',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
diff --git a/src/fft/ppm_fft_plan_3d_vec_bc2r_xy.f90 b/src/fft/ppm_fft_plan_3d_vec_bc2r_xy.f90
new file mode 100644
index 0000000..5655a6c
--- /dev/null
+++ b/src/fft/ppm_fft_plan_3d_vec_bc2r_xy.f90
@@ -0,0 +1,149 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_fft_plan_3d_vec_bc2r_xy
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! FFTW plan wrapper for 3d arrays, 2d complex to real
+ ! (backward) FFT in the xy directions
+ ! The routine does not work with fields that include ghost layers
+ !
+ ! For the C2R inverse transform FFTW always destroys the input data
+ !-------------------------------------------------------------------------
+#if __KIND == __SINGLE
+#define __ROUTINE ppm_fft_plan_3d_vec_bc2r_xy_s
+#define __PREC ppm_kind_single
+#elif __KIND == __DOUBLE
+#define __ROUTINE ppm_fft_plan_3d_vec_bc2r_xy_d
+#define __PREC ppm_kind_double
+#endif
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmplan,infield,outfield,info)
+ !!! FFTW plan wrapper for 3d arrays, 2d complex to real
+ !!! (backward) FFT in the xy directions
+ !!! The routine does not work with fields that include ghost layers
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_typedef
+ USE ppm_module_topo_get
+ USE ppm_module_write
+ USE ppm_module_data,ONLY:ppm_rank,ppm_kind_single,ppm_kind_double
+
+ IMPLICIT NONE
+
+ INCLUDE 'fftw3.f'
+
+ ! if debug check if dimensions are 2a 3b 5c 7d 11e 13f
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ !!!topology identifier of target
+ INTEGER,INTENT(IN) :: topoid
+ !!!id of the mesh
+ INTEGER,INTENT(IN) :: meshid
+ !!!ppm fft plan type
+ TYPE(ppm_fft_plan),INTENT(INOUT) :: ppmplan
+ !!!input field to fourier transform
+ !COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: infield
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: infield
+ !!!output field for the result of the fourier transform
+ !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: outfield
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: outfield
+ !!!Returns status, 0 upon success
+ INTEGER,INTENT(OUT) :: info
+ !in time perhaps an argument for alternate directions
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ REAL(__PREC) :: t0
+ INTEGER :: isub,isubl
+ INTEGER :: nsubs
+ INTEGER,DIMENSION(:),POINTER :: isublist
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_fft_plan',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_plan','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ nsubs = topology%nsublist
+ ALLOCATE(isublist(nsubs))
+ DO isub=1,nsubs
+ isublist(isub) = topology%isublist(isub)
+ ENDDO
+ mesh = topology%mesh(meshid)
+
+ !-------------------------------------------------------------------------
+ ! Setup parameters for this particular routine
+ !-------------------------------------------------------------------------
+ !the dimension of the FFT (1D/2D/3D)
+ ppmplan%rank=2
+ !the number of points along each direction of the piece to be transformed
+ ALLOCATE(ppmplan%nx(ppmplan%rank,nsubs))
+ !the direction of the transform
+ ppmplan%sign=FFTW_BACKWARD
+ !the method to setup the optimal plan
+ ppmplan%flag=FFTW_MEASURE
+ !the number of components to transform - 3 component vector
+ ppmplan%howmany=3
+ !the size of the input array - full size (assuming LBOUND=1 thus UBOUND)
+ ALLOCATE(ppmplan%inembed(ppmplan%rank))
+ ppmplan%inembed(1) = UBOUND(infield,2)
+ ppmplan%inembed(2) = UBOUND(infield,3)
+ !the size of the output array - full size (assuming LBOUND=1 thus UBOUND)
+ ALLOCATE(ppmplan%onembed(ppmplan%rank))
+ ppmplan%onembed(1) = UBOUND(outfield,2)
+ ppmplan%onembed(2) = UBOUND(outfield,3)
+ !istride tells how the same componenet data points are spaced in memory
+ !e.g. for 2/3 component vector istride = 2/3 or for scalar istride = 1
+ ppmplan%istride = 3
+ ppmplan%ostride = 3
+ !idist tells how multiple arrays are spaced in memory. I.e. a memory
+ !offset. e.g. vector components (idist=1) or scalar 2D arrays
+ !in 3D array(idist=NxNy)
+ ppmplan%idist = 1
+ ppmplan%odist = 1
+
+ !-------------------------------------------------------------------------
+ ! Allocate plan array
+ !-------------------------------------------------------------------------
+ IF(ASSOCIATED(ppmplan%plan)) THEN
+ DEALLOCATE(ppmplan%plan,stat=info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_plan','Failed to deallocate plan-array.',isub)
+ GOTO 9999
+ ENDIF
+ END IF
+ ALLOCATE(ppmplan%plan(nsubs))
+
+ DO isub=1,nsubs
+ isubl=isublist(isub)
+ !@ maybe the -1 needs to be removed when doing cell data
+ !we subtract the -1 to avoid the periodic vertex point
+ ppmplan%nx(1,isub) = mesh%nnodes(1,isubl)-1
+ ppmplan%nx(2,isub) = mesh%nnodes(2,isubl)-1
+
+ CALL dfftw_plan_many_dft_c2r(ppmplan%plan(isub),ppmplan%rank,&
+ & ppmplan%nx,ppmplan%howmany,infield(1,1,1,1,isub),ppmplan%inembed(1),&
+ & ppmplan%istride,ppmplan%idist,outfield(1,1,1,1,isub),&
+ & ppmplan%onembed(1),ppmplan%ostride,ppmplan%odist,ppmplan%flag)
+ END DO
+
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_fft_plan',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
diff --git a/src/fft/ppm_fft_plan_3d_vec_fc2c_z.f90 b/src/fft/ppm_fft_plan_3d_vec_fc2c_z.f90
new file mode 100644
index 0000000..45be5da
--- /dev/null
+++ b/src/fft/ppm_fft_plan_3d_vec_fc2c_z.f90
@@ -0,0 +1,145 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_fft_plan_3d_vec_fc2c_z
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! FFTW plan wrapper for 3d arrays, 1d complex to complex
+ ! (forward) FFT in the xy directions
+ ! The routine does not work with fields that include ghost layers
+ !-------------------------------------------------------------------------
+#if __KIND == __SINGLE
+#define __ROUTINE ppm_fft_plan_3d_vec_fc2c_z_s
+#define __PREC ppm_kind_single
+#elif __KIND == __DOUBLE
+#define __ROUTINE ppm_fft_plan_3d_vec_fc2c_z_d
+#define __PREC ppm_kind_double
+#endif
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmplan,infield,outfield,info)
+ !!! FFTW plan wrapper for 3d arrays, 1d complex to complex
+ !!! (forward) FFT in the xy directions
+ !!! The routine does not work with fields that include ghost layers
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_typedef
+ USE ppm_module_topo_get
+ USE ppm_module_write
+ USE ppm_module_data,ONLY:ppm_rank,ppm_kind_single,ppm_kind_double
+
+ IMPLICIT NONE
+
+ INCLUDE 'fftw3.f'
+
+ ! if debug check if dimensions are 2a 3b 5c 7d 11e 13f
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ !!!topology identifier of target
+ INTEGER,INTENT(IN) :: topoid
+ !!!id of the mesh
+ INTEGER,INTENT(IN) :: meshid
+ !!!ppm fft plan type
+ TYPE(ppm_fft_plan),INTENT(INOUT) :: ppmplan
+ !!!input field to fourier transform
+ !COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: infield
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: infield
+ !!!output field for the result of the fourier transform
+ !COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: outfield
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: outfield
+ !!!Returns status, 0 upon success
+ INTEGER,INTENT(OUT) :: info
+ !in time perhaps an argument for alternate directions
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ REAL(__PREC) :: t0
+ INTEGER :: isub,isubl
+ INTEGER :: nsubs
+ INTEGER,DIMENSION(:),POINTER :: isublist
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_fft_plan',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_plan','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ nsubs = topology%nsublist
+ ALLOCATE(isublist(nsubs))
+ DO isub=1,nsubs
+ isublist(isub) = topology%isublist(isub)
+ ENDDO
+ mesh = topology%mesh(meshid)
+
+ !-------------------------------------------------------------------------
+ ! Setup parameters for this particular routine
+ !-------------------------------------------------------------------------
+ !the dimension of the FFT (1D/2D/3D)
+ ppmplan%rank=1
+ !the number of points along each direction of the piece to be transformed
+ ALLOCATE(ppmplan%nx(ppmplan%rank,nsubs))
+ !the direction of the transform
+ ppmplan%sign=FFTW_FORWARD
+ !the method to setup the optimal plan
+ ppmplan%flag=FFTW_MEASURE
+ !the number of components to transform - 3 component vector
+ ppmplan%howmany=3
+ !the size of the input array - full size (assuming LBOUND=1 thus UBOUND)
+ ALLOCATE(ppmplan%inembed(ppmplan%rank))
+ ppmplan%inembed(1) = UBOUND(infield,4)
+ !the size of the output array - full size (assuming LBOUND=1 thus UBOUND)
+ ALLOCATE(ppmplan%onembed(ppmplan%rank))
+ ppmplan%onembed(1) = UBOUND(outfield,4)
+ !istride tells how the same componenet data points are spaced in memory
+ !e.g. for 2/3 component vector istride = 2/3 or for scalar istride = 1
+ ppmplan%istride = UBOUND(infield,2) *UBOUND(infield,3)*3
+ ppmplan%ostride = UBOUND(outfield,2)*UBOUND(outfield,3)*3
+ !idist tells how multiple arrays are spaced in memory. I.e. a memory
+ !offset. e.g. vector components (idist=1) or scalar 2D arrays in
+ !3D array(idist=NxNy)
+ ppmplan%idist = 1
+ ppmplan%odist = 1
+
+ !-------------------------------------------------------------------------
+ ! Allocate plan array
+ !-------------------------------------------------------------------------
+ IF(ASSOCIATED(ppmplan%plan)) THEN
+ DEALLOCATE(ppmplan%plan,stat=info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_plan','Failed to deallocate plan-array.',isub)
+ GOTO 9999
+ ENDIF
+ END IF
+ ALLOCATE(ppmplan%plan(nsubs))
+
+ DO isub=1,nsubs
+ isubl=isublist(isub)
+ !@ maybe the -1 needs to be removed when doing cell data
+ !we subtract the -1 to avoid the periodic vertex point
+ ppmplan%nx(1,isub) = mesh%nnodes(3,isubl)-1
+
+ CALL dfftw_plan_many_dft(ppmplan%plan(isub),ppmplan%rank,&
+ & ppmplan%nx(:,isub),ppmplan%howmany,infield(1,1,1,1,isub),&
+ & ppmplan%inembed(1),ppmplan%istride,ppmplan%idist,&
+ & outfield(1,1,1,1,isub),ppmplan%onembed(1),ppmplan%ostride,&
+ & ppmplan%odist,ppmplan%sign,ppmplan%flag)
+ END DO
+
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_fft_plan',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
diff --git a/src/fft/ppm_fft_plan_3d_vec_fr2c_xy.f90 b/src/fft/ppm_fft_plan_3d_vec_fr2c_xy.f90
new file mode 100644
index 0000000..79f4a81
--- /dev/null
+++ b/src/fft/ppm_fft_plan_3d_vec_fr2c_xy.f90
@@ -0,0 +1,148 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_fft_plan_3d_vec_fr2c_xy
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! FFTW plan wrapper for 3d arrays, 2d real to complex
+ ! (forward) FFT in the xy directions
+ ! The routine does not work with fields that include ghost layers
+ !-------------------------------------------------------------------------
+#if __KIND == __SINGLE
+#define __ROUTINE ppm_fft_plan_3d_vec_fr2c_xy_s
+#define __PREC ppm_kind_single
+#elif __KIND == __DOUBLE
+#define __ROUTINE ppm_fft_plan_3d_vec_fr2c_xy_d
+#define __PREC ppm_kind_double
+#endif
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmplan,infield,outfield,info)
+ !!! FFTW plan wrapper for 3d arrays, 2d real to complex
+ !!! (forward) FFT in the xy directions
+ !!! The routine does not work with fields that include ghost layers
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_typedef
+ USE ppm_module_topo_get
+ USE ppm_module_write
+ USE ppm_module_data,ONLY:ppm_rank,ppm_kind_single,ppm_kind_double
+
+ IMPLICIT NONE
+
+ INCLUDE 'fftw3.f'
+
+ ! if debug check if dimensions are 2a 3b 5c 7d 11e 13f
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ !!!topology identifier of target
+ INTEGER,INTENT(IN) :: topoid
+ !!!id of the mesh
+ INTEGER,INTENT(IN) :: meshid
+ !!!ppm fft plan type
+ TYPE(ppm_fft_plan),INTENT(INOUT) :: ppmplan
+ !!!input field to fourier transform
+ !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: infield
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: infield
+ !!!output field for the result of the fourier transform
+ !COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: outfield
+ COMPLEX(__PREC),DIMENSION(:,:,:,:,:),POINTER :: outfield
+ !!!Returns status, 0 upon success
+ INTEGER,INTENT(OUT) :: info
+ !in time perhaps an argument for alternate directions
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ REAL(__PREC) :: t0
+ INTEGER :: isub,isubl
+ INTEGER :: nsubs
+ INTEGER,DIMENSION(:),POINTER :: isublist
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_fft_plan',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_plan','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ nsubs = topology%nsublist
+ ALLOCATE(isublist(nsubs))
+ DO isub=1,nsubs
+ isublist(isub) = topology%isublist(isub)
+ ENDDO
+ mesh = topology%mesh(meshid)
+
+ !-------------------------------------------------------------------------
+ ! Setup parameters for this particular routine
+ !-------------------------------------------------------------------------
+ !the dimension of the FFT (1D/2D/3D)
+ ppmplan%rank=2
+ !the number of points along each direction of the piece to be transformed
+ ALLOCATE(ppmplan%nx(ppmplan%rank,nsubs))
+ !the direction of the transform
+ ppmplan%sign=FFTW_FORWARD
+ !the method to setup the optimal plan
+ ppmplan%flag=FFTW_MEASURE
+ !the number of components to transform - 3 component vector
+ ppmplan%howmany=3
+ !the size of the input array - full size (assuming LBOUND=1 thus UBOUND)
+ ALLOCATE(ppmplan%inembed(ppmplan%rank))
+ ppmplan%inembed(1) = UBOUND(infield,2)
+ ppmplan%inembed(2) = UBOUND(infield,3)
+ !the size of the output array - full size (assuming LBOUND=1 thus UBOUND)
+ ALLOCATE(ppmplan%onembed(ppmplan%rank))
+ ppmplan%onembed(1) = UBOUND(outfield,2)
+ ppmplan%onembed(2) = UBOUND(outfield,3)
+ !istride tells how the same componenet data points are spaced in memory
+ !e.g. for 2/3 component vector istride = 2/3 or for scalar istride = 1
+ ppmplan%istride = 3
+ ppmplan%ostride = 3
+ !idist tells how multiple arrays are spaced in memory. I.e. a memory
+ !offset. e.g. vector components (idist=1) or scalar 2D arrays in
+ !3D array(idist=NxNy)
+ ppmplan%idist = 1
+ ppmplan%odist = 1
+
+ !-------------------------------------------------------------------------
+ ! Allocate plan array
+ !-------------------------------------------------------------------------
+ IF(ASSOCIATED(ppmplan%plan)) THEN
+ DEALLOCATE(ppmplan%plan,stat=info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_fft_plan','Failed to deallocate plan-array.',isub)
+ GOTO 9999
+ ENDIF
+ END IF
+ ALLOCATE(ppmplan%plan(nsubs))
+
+ DO isub=1,nsubs
+ isubl=isublist(isub)
+ !@ maybe the -1 needs to be removed when doing cell data
+ !we subtract the -1 to avoid the periodic vertex point
+ ppmplan%nx(1,isub) = mesh%nnodes(1,isubl)-1
+ ppmplan%nx(2,isub) = mesh%nnodes(2,isubl)-1
+
+ CALL dfftw_plan_many_dft_r2c(ppmplan%plan(isub),ppmplan%rank,&
+ & ppmplan%nx(:,isub),ppmplan%howmany,infield(1,1,1,1,isub),&
+ & ppmplan%inembed(1),ppmplan%istride,ppmplan%idist,&
+ & outfield(1,1,1,1,isub),ppmplan%onembed(1),ppmplan%ostride,&
+ & ppmplan%odist,ppmplan%flag)
+ END DO
+
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_fft_plan',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
diff --git a/src/poisson/ppm_poisson_fd.f90 b/src/poisson/ppm_poisson_fd.f90
new file mode 100644
index 0000000..16a5cfd
--- /dev/null
+++ b/src/poisson/ppm_poisson_fd.f90
@@ -0,0 +1,144 @@
+ !-------------------------------------------------------------------------
+ ! ppm_poisson_fd.f90
+ !-------------------------------------------------------------------------
+ !-------------------------------------------------------------------------
+ #define __ROUTINE ppm_poisson_fd
+ #define __DIM 3
+ #define __NCOM 3
+ #define __ZEROSI (/0,0,0/)
+ SUBROUTINE __ROUTINE(topoid,meshid,fieldin,fieldout,dtype,info)
+
+ USE ppm_module_topo_get
+
+ IMPLICIT NONE
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ INTEGER, INTENT(IN) :: topoid
+ INTEGER, INTENT(IN) :: meshid
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldin
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldout
+ !INTEGER,DIMENSION(__DIM),INTENT(IN) :: gstw
+ INTEGER, INTENT(IN) :: dtype
+ INTEGER, INTENT(OUT) :: info
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ INTEGER,PARAMETER :: MK = __PREC
+ REAL(__PREC) :: t0
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+ REAL(__PREC) :: dx,dy,dz
+ REAL(__PREC) :: facx,facy,facz
+ INTEGER :: isub,isubl
+ INTEGER :: i,j,k
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_poisson_fd',t0,info)
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_init_predef','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ mesh = topology%mesh(meshid)
+
+ dx = (topology%max_physd(1)-topology%min_physd(1))/REAL(mesh%nm(1)-1)
+ dy = (topology%max_physd(2)-topology%min_physd(2))/REAL(mesh%nm(2)-1)
+ dz = (topology%max_physd(3)-topology%min_physd(3))/REAL(mesh%nm(3)-1)
+
+ !-----------------------------------------------------------------------
+ ! Do the finite difference calculation
+ !-----------------------------------------------------------------------
+ !-----------------------------------------------------------------------
+ ! Curl, 2nd order FD
+ !-----------------------------------------------------------------------
+ IF (dtype .EQ. ppm_poisson_drv_curl_fd2) THEN
+ facx = 1.0_MK/(2.0_MK*dx)
+ facy = 1.0_MK/(2.0_MK*dy)
+ facz = 1.0_MK/(2.0_MK*dz)
+
+ DO isub=1,topology%nsublist
+ isubl=topology%isublist(isub)
+ DO k=1,mesh%nnodes(3,isubl)
+ DO j=1,mesh%nnodes(2,isubl)
+ DO i=1,mesh%nnodes(1,isubl)
+ fieldout(1,i,j,k,isub) = &
+ & facz*(fieldin(2,i ,j ,k+1,isub)- &
+ & fieldin(2,i ,j ,k-1,isub)) &
+ & -facy*(fieldin(3,i ,j+1,k ,isub)- &
+ & fieldin(3,i ,j-1,k ,isub))
+ fieldout(2,i,j,k,isub) = &
+ & facx*(fieldin(3,i+1,j ,k ,isub)- &
+ & fieldin(3,i-1,j ,k ,isub)) &
+ & -facz*(fieldin(1,i ,j ,k+1,isub)- &
+ & fieldin(1,i ,j ,k-1,isub))
+ fieldout(3,i,j,k,isub) = &
+ & facy*(fieldin(1,i ,j+1,k ,isub)- &
+ & fieldin(1,i ,j-1,k ,isub)) &
+ & -facx*(fieldin(2,i+1,j ,k ,isub)- &
+ & fieldin(2,i-1,j ,k ,isub))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ !-----------------------------------------------------------------------
+ ! Curl, 4th order FD - untested
+ !-----------------------------------------------------------------------
+ ELSE IF (dtype .EQ. ppm_poisson_drv_curl_fd4) THEN
+ facx = 1.0_MK/(12.0_MK*dx)
+ facy = 1.0_MK/(12.0_MK*dy)
+ facz = 1.0_MK/(12.0_MK*dz)
+
+ DO isub=1,topology%nsublist
+ isubl=topology%isublist(isub)
+ DO k=1,mesh%nnodes(3,isubl)
+ DO j=1,mesh%nnodes(2,isubl)
+ DO i=1,mesh%nnodes(1,isubl)
+ fieldout(1,i,j,k,isub) = &
+ & facz*( -fieldin(2,i ,j ,k+2,isub) &
+ & +8.0_MK*fieldin(2,i ,j ,k+1,isub) &
+ & -8.0_MK*fieldin(2,i ,j ,k-1,isub) &
+ & +fieldin(2,i ,j ,k-2,isub)) &
+ & -facy*( -fieldin(3,i ,j+2,k ,isub) &
+ & +8.0_MK*fieldin(3,i ,j+1,k ,isub) &
+ & -8.0_MK*fieldin(3,i ,j-1,k ,isub) &
+ & +fieldin(3,i ,j-2,k ,isub))
+ fieldout(2,i,j,k,isub) = &
+ & facx*( -fieldin(3,i+2,j ,k ,isub) &
+ & +8.0_MK*fieldin(3,i+1,j ,k ,isub) &
+ & -8.0_MK*fieldin(3,i-1,j ,k ,isub) &
+ & +fieldin(3,i-2,j ,k ,isub)) &
+ & -facz*( -fieldin(1,i ,j ,k+2,isub) &
+ & +8.0_MK*fieldin(1,i ,j ,k+1,isub) &
+ & -8.0_MK*fieldin(1,i ,j ,k-1,isub) &
+ & +fieldin(1,i ,j ,k-2,isub))
+ fieldout(3,i,j,k,isub) = &
+ & facy*( -fieldin(1,i ,j+2,k ,isub) &
+ & +8.0_MK*fieldin(1,i ,j+1,k ,isub) &
+ & -8.0_MK*fieldin(1,i ,j-1,k ,isub) &
+ & +fieldin(1,i ,j-2,k ,isub)) &
+ & -facx*( -fieldin(2,i+2,j ,k ,isub) &
+ & +8.0_MK*fieldin(2,i+1,j ,k ,isub) &
+ & -8.0_MK*fieldin(2,i-1,j ,k ,isub) &
+ & +fieldin(2,i-2,j ,k ,isub))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_poisson_fd',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
+
diff --git a/src/poisson/ppm_poisson_init_predef.f90 b/src/poisson/ppm_poisson_init_predef.f90
new file mode 100644
index 0000000..0e3f695
--- /dev/null
+++ b/src/poisson/ppm_poisson_init_predef.f90
@@ -0,0 +1,414 @@
+ !-------------------------------------------------------------------------
+ ! ppm_poisson_init_predef.f90
+ !-------------------------------------------------------------------------
+ ! Routine to initialise the convolution using build in Greens function
+ !-------------------------------------------------------------------------
+#define __PREC ppm_kind_double
+!!#define __CMPLXDEF DCMPLX
+#define __DIM 3
+#define __ZEROSI (/0,0,0/)
+#define __ROUTINE ppm_poisson_init_predef
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmpoisson,fieldin,fieldout,green,info&
+ &,bc,derive)
+
+ USE ppm_module_mktopo
+ USE ppm_module_topo_get
+ USE ppm_module_mesh_define
+
+ IMPLICIT NONE
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ INTEGER, INTENT(IN) :: topoid
+ INTEGER, INTENT(IN) :: meshid
+ TYPE(ppm_poisson_plan),INTENT(INOUT) :: ppmpoisson
+ !@ strictly speaking fieldin is not being used in the init routine
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldin
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldout
+ !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: green
+ INTEGER, INTENT(IN) :: green
+ !!!flag to select build-in Greens functions:
+ !!!ppm_poisson_grn_pois_per - Poisson equation, periodic boundaries
+ !!!ppm_poisson_grn_pois_fre - Poisson equation, freespace boundaries (not implemented)
+ !!!ppm_poisson_grn_reprojec - Do vorticity reprojection to kill divergence
+ INTEGER, INTENT(OUT) :: info
+ INTEGER,INTENT(IN),OPTIONAL :: bc
+ !!!boundary condition for the convolution. Can be on of the following:
+ !!!ppm_poisson_grn_pois_per, ppm_poisson_grn_pois_fre.
+ !!!One could argue that this is redundant in the build-in case
+ !!!@ Isnt this redundant because of green?
+ INTEGER,INTENT(IN),OPTIONAL :: derive
+ !!!flag to toggle various derivatives of the solution (not to be used with
+ !!!green=ppm_poisson_grn_reprojec):
+ !!!ppm_poisson_drv_none
+ !!!ppm_poisson_drv_curl_fd (not implemented)
+ !!!ppm_poisson_drv_grad_fd (not implemented)
+ !!!ppm_poisson_drv_lapl_fd (not implemented)
+ !!!ppm_poisson_drv_curl_sp (not implemented)
+ !!!ppm_poisson_drv_grad_sp (not implemented)
+ !!!ppm_poisson_drv_lapl_sp (not implemented)
+ !!!The spectral derivatives can only be computed in this routine. Since
+ !!!the flag exists finite difference derivatives have also been included.
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ REAL(__PREC) :: t0
+ REAL(__PREC),DIMENSION(:,:),POINTER :: xp !particle positions
+ TYPE(ppm_t_topo),POINTER :: topology,topologyxy,topologyxyc,topologyz
+ TYPE(ppm_t_equi_mesh) :: mesh!!,meshtmp1,meshtmp2 !@
+ INTEGER ,DIMENSION(__DIM) :: indl,indu
+ INTEGER,PARAMETER :: MK = __PREC
+ REAL(__PREC),PARAMETER :: PI=ACOS(-1.0_MK) !@ use ppm pi
+ !factor for the Green's function, including FFT normalization
+ REAL(__PREC) :: normfac
+ INTEGER :: i,j,k
+ INTEGER :: kx,ky,kz
+ INTEGER :: isubl,isub
+ INTEGER,DIMENSION(__DIM*2) :: bcdef
+ INTEGER :: assigning
+ INTEGER :: decomposition
+ INTEGER,SAVE :: ttopoid
+ INTEGER :: tmeshid
+
+ REAL(__PREC),DIMENSION(__DIM) :: tmpmin,tmpmax
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_poisson_init_predef',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Investigate optional arguments, setup routine accordingly
+ ! !@ Also check if the input/output and derivatives match
+ !-------------------------------------------------------------------------
+ IF (green .EQ. ppm_poisson_grn_pois_per) THEN
+ ppmpoisson%case = ppm_poisson_grn_pois_per
+ ELSE IF (green .EQ. ppm_poisson_grn_pois_fre) THEN
+ ppmpoisson%case = ppm_poisson_grn_pois_fre
+ ELSE IF (green .EQ. ppm_poisson_grn_reprojec) THEN
+ ppmpoisson%case = ppm_poisson_grn_reprojec
+ ENDIF
+
+ IF (PRESENT(derive)) THEN
+ !Create temporary arrays if necessary
+ IF (derive .EQ. ppm_poisson_drv_curl_fd2 .AND.&
+ & ppmpoisson%case .EQ. ppm_poisson_grn_pois_per) THEN
+ ppmpoisson%derivatives = derive
+ ALLOCATE(ppmpoisson%drv_vr(LBOUND(fieldout,1):UBOUND(fieldout,1),&
+ & LBOUND(fieldout,2):UBOUND(fieldout,2),&
+ & LBOUND(fieldout,3):UBOUND(fieldout,3),&
+ & LBOUND(fieldout,4):UBOUND(fieldout,4),&
+ & LBOUND(fieldout,5):UBOUND(fieldout,5)))
+ ELSE IF (derive .EQ. ppm_poisson_drv_none) THEN
+ ppmpoisson%derivatives = ppm_poisson_drv_none
+ ELSE
+ CALL ppm_write(ppm_rank,'ppm_poisson_init_predef','Undefined derivation input.',isub)
+ GOTO 9999
+ ENDIF
+ ELSE
+ ppmpoisson%derivatives = ppm_poisson_drv_none
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Get topology and mesh values
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_init_predef','Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ !nsubs = topology%nsublist
+ mesh = topology%mesh(meshid)
+ !Copy size of global mesh
+ ppmpoisson%nmxy (1) = mesh%nm(1)
+ ppmpoisson%nmxy (2) = mesh%nm(2)
+ ppmpoisson%nmxy (3) = mesh%nm(3)
+ !ppmpoisson%nmxyc(1) = (mesh%nm(1)-1)/2+1
+ ppmpoisson%nmxyc(1) = (mesh%nm(1)) !tmp until meshid>1 works
+ ppmpoisson%nmxyc(2) = mesh%nm(2)
+ ppmpoisson%nmxyc(3) = mesh%nm(3)
+ !ppmpoisson%nmz (1) = (mesh%nm(1)-1)/2+1!/2+1 !trying to reduce the number of points in matrix+1
+ ppmpoisson%nmz (1) = (mesh%nm(1)) !tmp until meshid>1 works
+ ppmpoisson%nmz (2) = mesh%nm(2)
+ ppmpoisson%nmz (3) = mesh%nm(3)
+
+
+ !-------------------------------------------------------------------------
+ ! Create new slab topology !@sofar periodic
+ !-------------------------------------------------------------------------
+ ttopoid = 0
+ tmeshid = -1
+ decomposition = ppm_param_decomp_xy_slab
+ assigning = ppm_param_assign_internal
+ bcdef = ppm_param_bcdef_periodic
+ tmpmin = topology%min_physd
+ tmpmax = topology%max_physd
+ !!tmpmin2 = topology%min_physd
+ !!tmpmax2 = topology%max_physd
+
+ CALL ppm_mktopo(ttopoid,tmeshid,xp,0,&
+ & decomposition,assigning,&
+ !& topology%min_physd,topology%max_physd,bcdef,&
+ & tmpmin,tmpmax,bcdef,&
+ & __ZEROSI,ppmpoisson%costxy,ppmpoisson%istartxy,ppmpoisson%ndataxy,&
+ & ppmpoisson%nmxy,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_init_predef','Failed to create xy-topology.',isub)
+ GOTO 9999
+ ENDIF
+ ppmpoisson%topoidxy = ttopoid
+ ppmpoisson%meshidxy = tmeshid
+
+ !-------------------------------------------------------------------------
+ ! Create complex slab mesh
+ ! !@ not used until meshid>1 works
+ !-------------------------------------------------------------------------
+ ttopoid = ppmpoisson%topoidxy
+ tmeshid = -1
+ CALL ppm_mesh_define(ttopoid,tmeshid,&
+ & ppmpoisson%nmxy,ppmpoisson%istartxyc,ppmpoisson%ndataxyc,info) !@nmxyC
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_init_predef','Failed to create complex xy mesh definition.',isub)
+ GOTO 9999
+ ENDIF
+ ppmpoisson%meshidxyc = tmeshid
+ !write(*,*) ppmpoisson%meshidxy,tmeshid,ttopoid,ppmpoisson%meshidxyc
+ !!write(*,*) ppmpoisson%istartxyc,ppmpoisson%ndataxyc
+
+
+ !-------------------------------------------------------------------------
+ ! Create new pencil topology !@sofar periodic
+ !-------------------------------------------------------------------------
+ ttopoid = 0
+ tmeshid = -1
+ bcdef = ppm_param_bcdef_periodic
+ assigning = ppm_param_assign_internal
+ decomposition = ppm_param_decomp_zpencil
+ !tmpmin = topology%min_physd
+ !tmpmax = topology%max_physd
+
+ CALL ppm_mktopo(ttopoid,tmeshid,xp,0,&
+ & decomposition,assigning,&
+ !& topology%min_physd,topology%max_physd,bcdef,&
+ & tmpmin,tmpmax,bcdef,&
+ & __ZEROSI,ppmpoisson%costz,ppmpoisson%istartz,ppmpoisson%ndataz,&
+ & ppmpoisson%nmz,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_init_predef','Failed to create z-topology.',isub)
+ GOTO 9999
+ ENDIF
+ ppmpoisson%topoidz = ttopoid
+ ppmpoisson%meshidz = tmeshid
+
+ !-------------------------------------------------------------------------
+ ! Get the new topologies (slabs and pencils)
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(ppmpoisson%topoidxy,topologyxy,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_init_predef','Failed to get xy topology.',isub)
+ GOTO 9999
+ ENDIF
+
+ !!CALL ppm_topo_get(ppmpoisson%topoidxyc,topologyxyc,info)
+ !!IF (info .NE. 0) THEN
+ !!CALL ppm_write(ppm_rank,'ppm_poisson_init_predef','Failed to get complex xy topology.',isub)
+ !!GOTO 9999
+ !!ENDIF
+
+ CALL ppm_topo_get(ppmpoisson%topoidz,topologyz,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_init_predef','Failed to get z topology.',isub)
+ GOTO 9999
+ ENDIF
+
+
+ !-------------------------------------------------------------------------
+ ! Copy data to ppm_poisson type
+ !-------------------------------------------------------------------------
+ ppmpoisson%nsublistxy = topologyxy %nsublist
+ !!ppmpoisson%nsublistxyc = topologyxyc%nsublist
+ ppmpoisson%nsublistz = topologyz %nsublist
+
+ ALLOCATE(ppmpoisson%isublistxy (ppmpoisson%nsublistxy))
+ !!ALLOCATE(ppmpoisson%isublistxyc(ppmpoisson%nsublistxyc))
+ ALLOCATE(ppmpoisson%isublistz (ppmpoisson%nsublistz))
+
+ DO isub=1,ppmpoisson%nsublistxy
+ ppmpoisson%isublistxy(isub) = topologyxy%isublist(isub)
+ ENDDO
+ !!DO isub=1,ppmpoisson%nsublistxyc
+ !!ppmpoisson%isublistxyc(isub) = topologyxyc%isublist(isub)
+ !!ENDDO
+ DO isub=1,ppmpoisson%nsublistz
+ ppmpoisson%isublistz(isub) = topologyz%isublist(isub)
+ ENDDO
+
+ !-------------------------------------------------------------------------
+ ! Set and get minimum and maximum indicies of xy slabs
+ !-------------------------------------------------------------------------
+ indl(1) = 1
+ indl(2) = 1
+ indl(3) = 1
+ indu(1) = 0
+ indu(2) = 0
+ indu(3) = 0
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl = ppmpoisson%isublistxy(isub)
+ indu(1) = MAX(indu(1),ppmpoisson%ndataxy(1,isubl))
+ indu(2) = MAX(indu(2),ppmpoisson%ndataxy(2,isubl))
+ indu(3) = MAX(indu(3),ppmpoisson%ndataxy(3,isubl))
+ ENDDO
+
+ !-------------------------------------------------------------------------
+ ! Allocate real xy slabs
+ !-------------------------------------------------------------------------
+ ALLOCATE(ppmpoisson%fldxyr(__DIM,&
+ & indl(1):indu(1),indl(2):indu(2),indl(3):indu(3),&
+ & 1:ppmpoisson%nsublistxy),stat=info)
+
+ !!ALLOCATE(ppmpoisson%fldxyc(__DIM,&
+ !!& indl(1):indu(1),indl(2):indu(2),indl(3):indu(3),&
+ !!& 1:ppmpoisson%nsublistxy),stat=info)
+
+ !-------------------------------------------------------------------------
+ ! Set and get minimum and maximum indicies of COMPLEX xy slabs
+ !-------------------------------------------------------------------------
+ indl(1) = 1
+ indl(2) = 1
+ indl(3) = 1
+ indu(1) = 0
+ indu(2) = 0
+ indu(3) = 0
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl = ppmpoisson%isublistxy(isub)
+ indu(1) = MAX(indu(1),ppmpoisson%ndataxyc(1,isubl))
+ indu(2) = MAX(indu(2),ppmpoisson%ndataxyc(2,isubl))
+ indu(3) = MAX(indu(3),ppmpoisson%ndataxyc(3,isubl))
+ ENDDO
+
+ !-------------------------------------------------------------------------
+ ! Allocate real and complex xy slabs
+ !-------------------------------------------------------------------------
+ !!ALLOCATE(ppmpoisson%fldxyr(__DIM,&
+ !!& indl(1):indu(1),indl(2):indu(2),indl(3):indu(3),&
+ !!& 1:ppmpoisson%nsublistxy),stat=info)
+!!
+ ALLOCATE(ppmpoisson%fldxyc(__DIM,&
+ & indl(1):indu(1),indl(2):indu(2),indl(3):indu(3),&
+ & 1:ppmpoisson%nsublistxy),stat=info)
+
+ !-------------------------------------------------------------------------
+ ! Set and get minimum and maximum indicies of z pencils
+ !-------------------------------------------------------------------------
+ indl(1) = 1
+ indl(2) = 1
+ indl(3) = 1
+ indu(1) = 0
+ indu(2) = 0
+ indu(3) = 0
+ DO isub=1,ppmpoisson%nsublistz
+ isubl = ppmpoisson%isublistz(isub)
+ indu(1) = MAX(indu(1),ppmpoisson%ndataz(1,isubl))
+ indu(2) = MAX(indu(2),ppmpoisson%ndataz(2,isubl))
+ indu(3) = MAX(indu(3),ppmpoisson%ndataz(3,isubl))
+ ENDDO
+
+ !-------------------------------------------------------------------------
+ ! Allocate two complex z pencils + Greens fcn array !@check return vars.
+ !-------------------------------------------------------------------------
+ ALLOCATE(ppmpoisson%fldzc1(__DIM,&
+ & indl(1):indu(1),indl(2):indu(2),indl(3):indu(3),&
+ & 1:ppmpoisson%nsublistz),stat=info)
+
+ ALLOCATE(ppmpoisson%fldzc2(__DIM,&
+ & indl(1):indu(1),indl(2):indu(2),indl(3):indu(3),&
+ & 1:ppmpoisson%nsublistz),stat=info)
+
+ IF (green .EQ. ppm_poisson_grn_pois_per) THEN
+ ALLOCATE(ppmpoisson%fldgrnr(&
+ & indl(1):indu(1),indl(2):indu(2),indl(3):indu(3),&
+ & 1:ppmpoisson%nsublistz),stat=info)
+ ELSE IF (green .EQ. ppm_poisson_grn_pois_fre) THEN
+ !@ freespace solutions are not even remotely implemented
+ ALLOCATE(ppmpoisson%fldgrnc(3,&
+ & indl(1):indu(1),indl(2):indu(2),indl(3):indu(3),&
+ & 1:ppmpoisson%nsublistz),stat=info)
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Set up xy FFT plans
+ !-------------------------------------------------------------------------
+ CALL ppm_fft_forward_2d(ppmpoisson%topoidxy,ppmpoisson%meshidxy,&
+ & ppmpoisson%planfxy,ppmpoisson%fldxyr,&
+ & ppmpoisson%fldxyc,info)
+
+ CALL ppm_fft_backward_2d(ppmpoisson%topoidxy,ppmpoisson%meshidxy,&
+ & ppmpoisson%planbxy,ppmpoisson%fldxyc,&
+ & ppmpoisson%fldxyr,info)
+
+ !-------------------------------------------------------------------------
+ ! Set up z FFT plans
+ !-------------------------------------------------------------------------
+ CALL ppm_fft_forward_1d(ppmpoisson%topoidz,ppmpoisson%meshidz,&
+ & ppmpoisson%planfz,ppmpoisson%fldzc1,&
+ & ppmpoisson%fldzc2,info)
+
+ CALL ppm_fft_backward_1d(ppmpoisson%topoidz,ppmpoisson%meshidz,&
+ & ppmpoisson%planbz,ppmpoisson%fldzc2,&
+ & ppmpoisson%fldzc1,info)
+
+ !-------------------------------------------------------------------------
+ ! Compute Greens function. Analytic, periodic
+ !
+ ! (d2_/dx2 + d2_/dy2 + d2_/dz2)psi = -omega =>
+ ! -4*pi2(kx2 + ky2 + kz2)PSI = -OMEGA =>
+ ! PSI = 1/(4*pi2)*1/(kx2 + ky2 + kz2)OMEGA
+ !-------------------------------------------------------------------------
+ ! Scaling the spectral koefficients... and normalisation of FFTs
+ normfac = -1.0_MK/(4.0_MK*PI*PI * &
+ & REAL((mesh%nm(1)-1)*(mesh%nm(2)-1)*(mesh%nm(3)-1),MK))
+ IF (green .EQ. ppm_poisson_grn_pois_per) THEN
+ DO isub=1,ppmpoisson%nsublistz
+ isubl=ppmpoisson%isublistz(isub)
+ DO k=1,ppmpoisson%ndataz(3,isubl)
+ DO j=1,ppmpoisson%ndataz(2,isubl)
+ DO i=1,ppmpoisson%ndataz(1,isubl)
+ kx = i-1 + (ppmpoisson%istartz(1,isubl)-1)
+ ky = j-1 + (ppmpoisson%istartz(2,isubl)-1)
+ kz = k-1 + (ppmpoisson%istartz(3,isubl)-1)
+ !The integer division depends on mesh%nm to include N+1 points:
+ !This is a nasty way to do this but it is only done once so...:
+ IF (kx .GT. (mesh%nm(1)-1)/2) kx = kx-(mesh%nm(1)-1)
+ IF (ky .GT. (mesh%nm(2)-1)/2) ky = ky-(mesh%nm(2)-1)
+ IF (kz .GT. (mesh%nm(3)-1)/2) kz = kz-(mesh%nm(3)-1)
+ ppmpoisson%fldgrnr(i,j,k,isub) = &
+ & normfac/(REAL(kx*kx+ky*ky+kz*kz,__PREC))
+ !Take care of singularity
+ !This is nasty as well
+ IF ((kx*kx+ky*ky+kz*kz) .EQ. 0) THEN
+ ppmpoisson%fldgrnr(i,j,k,isub) = 0.0_MK
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ END IF
+
+ !-------------------------------------------------------------------------
+ ! Or alternatively FFT real Green from input
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Deallocate fields? !@
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_poisson_init_predef',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
+
diff --git a/src/poisson/ppm_poisson_solve.f90 b/src/poisson/ppm_poisson_solve.f90
new file mode 100644
index 0000000..6741625
--- /dev/null
+++ b/src/poisson/ppm_poisson_solve.f90
@@ -0,0 +1,712 @@
+ !-------------------------------------------------------------------------
+ ! ppm_poisson_solve.f90
+ !-------------------------------------------------------------------------
+ !@mapping calls can be cleaned up
+ !-------------------------------------------------------------------------
+ #define __ROUTINE ppm_poisson_solve
+ #define __DIM 3
+ #define __NCOM 3
+ #define __ZEROSI (/0,0,0/)
+ !#define __NOPE
+ SUBROUTINE __ROUTINE(topoid,meshid,ppmpoisson,fieldin,fieldout,gstw,info,&
+ & tmpcase)
+
+ USE ppm_module_map_field
+ USE ppm_module_map_field_global
+ USE ppm_module_map
+
+
+ IMPLICIT NONE
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ INTEGER, INTENT(IN) :: topoid
+ INTEGER, INTENT(IN) :: meshid
+ TYPE(ppm_poisson_plan),INTENT(INOUT) :: ppmpoisson
+ !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: fieldin
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldin
+ !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: fieldout
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldout
+ INTEGER,DIMENSION(__DIM),INTENT(IN) :: gstw
+ INTEGER, INTENT(OUT) :: info
+ INTEGER,OPTIONAL,INTENT(IN) :: tmpcase
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ INTEGER,PARAMETER :: MK = __PREC
+ REAL(__PREC) :: t0
+ INTEGER :: isub,isubl
+ INTEGER :: i,j,k
+ INTEGER :: info2
+ INTEGER :: presentcase
+ REAL(__PREC) :: wdotk
+ INTEGER :: gi,gj,gk
+ REAL(__PREC) :: kx,ky,kz
+ REAL(__PREC) :: normfac
+#ifndef __NOPE
+INTEGER :: trank !@
+trank =0
+#endif
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_poisson_solve',t0,info)
+
+ !-------------------------------------------------------------------------
+ ! Check if we run a different/temporary case
+ !-------------------------------------------------------------------------
+ IF (PRESENT(tmpcase)) THEN
+ presentcase = tmpcase
+ ELSE
+ presentcase = ppmpoisson%case
+ ENDIF
+
+
+ !-------------------------------------------------------------------------
+ ! Perhaps allocate (and deallocate) arrays !@
+ !-------------------------------------------------------------------------
+
+
+ !------------------------------------------------------------------------
+ ! Map data globally to the slabs (XY)
+ !------------------------------------------------------------------------
+ !Initialise
+ CALL ppm_map_field_global(&
+ & topoid, &
+ & ppmpoisson%topoidxy, &
+ & meshid, &
+ & ppmpoisson%meshidxy,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_solve','Failed to initialise field mapping.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Push the data
+ CALL ppm_map_field_push(&
+ & topoid, &
+ & meshid,fieldin,__NCOM,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to push vector field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Send
+ CALL ppm_map_field_send(info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to send field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Retrieve
+ CALL ppm_map_field_pop(&
+ & ppmpoisson%topoidxy, &
+ & ppmpoisson%meshidxy,ppmpoisson%fldxyr, &
+ & __NCOM,__ZEROSI,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !!ppmpoisson%fldxyc= 0.0_MK
+ !!ppmpoisson%fldzc1= 0.0_MK
+ !!ppmpoisson%fldzc2= 0.0_MK
+ !-----------------------------------------------------------------------
+ ! Do slab FFT (XY) - use the non-reduced topology
+ !-----------------------------------------------------------------------
+ CALL ppm_fft_execute_2d(ppmpoisson%topoidxy,&
+ & ppmpoisson%meshidxy, ppmpoisson%planfxy, &
+ & ppmpoisson%fldxyr, ppmpoisson%fldxyc, &
+ & info)
+
+ !DO isub=1,ppmpoisson%nsublistxy
+ !isubl=ppmpoisson%isublistxy(isub)
+ !!Copy x/y boundaries
+ !DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ !DO i=1,ppmpoisson%ndataxy(1,isubl)
+ !ppmpoisson%fldxyc(1,i,1,k,isub) = ppmpoisson%fldxyc(1,i,ppmpoisson%ndataxy(2,isubl),k,isub)
+ !END DO
+ !DO j=1,ppmpoisson%ndataxy(2,isubl)
+ !ppmpoisson%fldxyc(1,1,j,k,isub) = ppmpoisson%fldxyc(1,ppmpoisson%ndataxy(1,isubl),j,k,isub)
+ !END DO
+ !END DO
+ !END DO
+
+ !!ppmpoisson%fldxyr = 0.0_MK
+ !!fieldin = 0.0_MK
+ #ifdef __NOPE
+ if (ppm_rank .EQ. trank) THEN
+ write(*,*) 'Rfftx12', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)-1
+ write(*,'(A,E12.4,A,$)') ' (',ppmpoisson%fldxyr(1,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ write(*,*) 'Rffty12', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)-1
+ write(*,'(A,E12.4,A,$)') ' (',ppmpoisson%fldxyr(2,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ write(*,*) 'Rfftz12', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)-1
+ write(*,'(A,E12.4,A,$)') ' (',ppmpoisson%fldxyr(3,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ ENDIF
+ #endif
+
+
+ #ifdef __NOPE
+ if (ppm_rank .EQ. trank) THEN
+ write(*,*) 'fftx12', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)-1
+ write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldxyc(1,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ write(*,*) 'ffty12', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)-1
+ write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldxyc(2,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ write(*,*) 'fftz12', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)-1
+ write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldxyc(3,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ ENDIF
+ #endif
+
+ !#ifdef __NOPE
+ !-----------------------------------------------------------------------
+ ! Map to the pencils (Z)
+ !-----------------------------------------------------------------------
+ !Initialise
+ CALL ppm_map_field_global(&
+ & ppmpoisson%topoidxy, &
+ & ppmpoisson%topoidz, &
+ !& ppmpoisson%meshidxyc, & !@to be used when meshid>1 works
+ & ppmpoisson%meshidxy, &
+ & ppmpoisson%meshidz,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to initialise field mapping.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Push the data
+ CALL ppm_map_field_push(&
+ & ppmpoisson%topoidxy, &
+ & ppmpoisson%meshidxyc,ppmpoisson%fldxyc,__NCOM,info)
+ !& ppmpoisson%meshidxyc,ppmpoisson%fldxyc,__NCOM,info)!@to be used when meshid>1 works
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to push vector field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Send
+ CALL ppm_map_field_send(info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to send field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Retrieve
+ CALL ppm_map_field_pop(&
+ & ppmpoisson%topoidz, &
+ & ppmpoisson%meshidz,ppmpoisson%fldzc1, &
+ & __NCOM,__ZEROSI,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !-----------------------------------------------------------------------
+ ! Do pencil FFT (Z)
+ !-----------------------------------------------------------------------
+ CALL ppm_fft_execute_1d(ppmpoisson%topoidz,&
+ & ppmpoisson%meshidz, ppmpoisson%planfz, &
+ & ppmpoisson%fldzc1, ppmpoisson%fldzc2, &
+ & info)
+
+ #ifdef __NOPE
+ if (ppm_rank .EQ. trank) THEN
+ write(*,*)
+ write(*,*)
+ !ppmpoisson%fldzc1 = 0.0_MK
+ write(*,*) 'fftx123', ppm_rank
+ DO isub=1,ppmpoisson%nsublistz
+ isubl=ppmpoisson%isublistz(isub)
+ DO k=1,ppmpoisson%ndataz(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataz(1,isubl)-1
+ DO j=1,ppmpoisson%ndataz(2,isubl)-1
+ write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldzc2(1,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ write(*,*) 'ffty123', ppm_rank
+ DO isub=1,ppmpoisson%nsublistz
+ isubl=ppmpoisson%isublistz(isub)
+ DO k=1,ppmpoisson%ndataz(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataz(1,isubl)-1
+ DO j=1,ppmpoisson%ndataz(2,isubl)-1
+ write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldzc2(2,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ write(*,*) 'fftz123', ppm_rank
+ DO isub=1,ppmpoisson%nsublistz
+ isubl=ppmpoisson%isublistz(isub)
+ DO k=1,ppmpoisson%ndataz(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataz(1,isubl)-1
+ DO j=1,ppmpoisson%ndataz(2,isubl)-1
+ write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldzc2(3,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+
+ write(*,*)
+ write(*,*)
+ write(*,*) 'green', ppm_rank
+ DO isub=1,ppmpoisson%nsublistz
+ isubl=ppmpoisson%isublistz(isub)
+ DO k=1,ppmpoisson%ndataz(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataz(1,isubl)-1
+ DO j=1,ppmpoisson%ndataz(2,isubl)-1
+ write(*,'(E12.4,$)') REAL(ppmpoisson%fldgrnr(i,j,k,isub))
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ ENDIF
+ #endif
+
+ !-----------------------------------------------------------------------
+ ! Apply the Green's function
+ ! ... depending on the case
+ !-----------------------------------------------------------------------
+ !!kill?: ppmpoisson%fldgrnc = -10000000.0_MK
+ !!kill?: DO isub=1,ppmpoisson%nsublistz
+ !!kill?: isubl=ppmpoisson%isublistz(isub)
+ !!kill?: DO k=1,ppmpoisson%ndataz(3,isubl)
+ !!kill?: DO j=1,ppmpoisson%ndataz(2,isubl)
+ !!kill?: DO i=1,ppmpoisson%ndataz(1,isubl)
+ !!kill?: ppmpoisson%fldgrnc(1,i,j,k,isub) = ppmpoisson%fldzc2(1,i,j,k,isub)
+ !!kill?: ppmpoisson%fldgrnc(2,i,j,k,isub) = ppmpoisson%fldzc2(2,i,j,k,isub)
+ !!kill?: ppmpoisson%fldgrnc(3,i,j,k,isub) = ppmpoisson%fldzc2(3,i,j,k,isub)
+ !!kill?: !ppmpoisson%fldzc2(1,i,j,k,isub) = ppmpoisson%fldgrnr(i,j,k,isub)
+ !!kill?: !ppmpoisson%fldzc2(2,i,j,k,isub) = ppmpoisson%fldgrnr(i,j,k,isub)
+ !!kill?: !ppmpoisson%fldzc2(3,i,j,k,isub) = ppmpoisson%fldgrnr(i,j,k,isub)
+ !!kill?: !!ppmpoisson%fldzc2(1,i,j,k,isub) = CMPLX((REAL(i*j*k)),0.0)
+ !!kill?: !!ppmpoisson%fldzc2(2,i,j,k,isub) = CMPLX((REAL(i*j*k)),0.0)
+ !!kill?: !!ppmpoisson%fldzc2(3,i,j,k,isub) = CMPLX((REAL(i*j*k)),0.0)
+ !!kill?: !fieldout(1,i,j,k,isub) = REAL(ppmpoisson%fldgrnr(i,j,k,isub))
+ !!kill?: !fieldout(2,i,j,k,isub) = REAL(ppmpoisson%fldgrnr(i,j,k,isub))
+ !!kill?: !fieldout(3,i,j,k,isub) = REAL(ppmpoisson%fldgrnr(i,j,k,isub))
+ !!kill?: ENDDO
+ !!kill?: ENDDO
+ !!kill?: ENDDO
+ !!kill?: ENDDO
+
+ IF (presentcase .EQ. ppm_poisson_grn_pois_per) THEN
+ DO isub=1,ppmpoisson%nsublistz
+ isubl=ppmpoisson%isublistz(isub)
+ DO k=1,ppmpoisson%ndataz(3,isubl)
+ DO j=1,ppmpoisson%ndataz(2,isubl)
+ DO i=1,ppmpoisson%ndataz(1,isubl)
+ ppmpoisson%fldzc2(1,i,j,k,isub) = ppmpoisson%fldgrnr( i,j,k,isub)*&
+ & ppmpoisson%fldzc2(1,i,j,k,isub)
+ ppmpoisson%fldzc2(2,i,j,k,isub) = ppmpoisson%fldgrnr( i,j,k,isub)*&
+ & ppmpoisson%fldzc2(2,i,j,k,isub)
+ ppmpoisson%fldzc2(3,i,j,k,isub) = ppmpoisson%fldgrnr( i,j,k,isub)*&
+ & ppmpoisson%fldzc2(3,i,j,k,isub)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ELSE IF (presentcase .EQ. ppm_poisson_grn_reprojec) THEN
+ !remember to normalize the FFT
+ normfac = 1.0_MK/ REAL((ppmpoisson%nmz(1)-1)* &
+ & (ppmpoisson%nmz(2)-1)* &
+ & (ppmpoisson%nmz(3)-1),MK)
+ write(*,*) 'reprojection ppm style'
+ DO isub=1,ppmpoisson%nsublistz
+ isubl=ppmpoisson%isublistz(isub)
+ DO k=1,ppmpoisson%ndataz(3,isubl)
+ gk = k - 1 + (ppmpoisson%istartz(3,isubl)-1)
+ IF (gk .GT. (ppmpoisson%nmz(3)-1)/2) gk = gk-(ppmpoisson%nmz(3)-1)
+ kz = REAL(gk,MK)
+ DO j=1,ppmpoisson%ndataz(2,isubl)
+ gj = j - 1 + (ppmpoisson%istartz(2,isubl)-1)
+ IF (gj .GT. (ppmpoisson%nmz(2)-1)/2) gj = gj-(ppmpoisson%nmz(2)-1)
+ ky = REAL(gj,MK)
+ DO i=1,ppmpoisson%ndataz(1,isubl)
+ gi = i - 1 + (ppmpoisson%istartz(1,isubl)-1)
+ IF (gi .GT. (ppmpoisson%nmz(1)-1)/2) gi = gi-(ppmpoisson%nmz(1)-1)
+ kx = REAL(gi,MK)
+
+ IF (gi .EQ. 0 .AND. gj .EQ. 0 .AND. gk .EQ. 0) THEN
+ wdotk = 0.0_mk
+ ELSE
+ wdotk = (ppmpoisson%fldzc2(1,i,j,k,isub) * kx + &
+ & ppmpoisson%fldzc2(2,i,j,k,isub) * ky + &
+ & ppmpoisson%fldzc2(3,i,j,k,isub) * kz) / &
+ & (kx*kx+ky*ky+kz*kz)
+ ENDIF
+
+ !!write(*,*) ppmpoisson%fldzc2(:,i,j,k,isub),wdotk
+
+ ppmpoisson%fldzc2(1,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(1,i,j,k,isub) - wdotk*kx)*normfac
+ ppmpoisson%fldzc2(2,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(2,i,j,k,isub) - wdotk*ky)*normfac
+ ppmpoisson%fldzc2(3,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(3,i,j,k,isub) - wdotk*kz)*normfac
+ !!write(*,*) ppmpoisson%fldzc2(:,i,j,k,isub)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+
+ !write(*,*) ppmpoisson%fldzc2
+ !-----------------------------------------------------------------------
+ ! IFFT pencil (Z)
+ !-----------------------------------------------------------------------
+ CALL ppm_fft_execute_1d(ppmpoisson%topoidz,&
+ & ppmpoisson%meshidz, ppmpoisson%planbz, &
+ & ppmpoisson%fldzc2, ppmpoisson%fldzc1, &
+ & info)
+ !!CALL ppm_fft_normalize(ppmpoisson%topoidz,&
+ !!& ppmpoisson%meshidz, ppmpoisson%planbz,__ZEROSI,&
+ !!& ppmpoisson%fldzc1, &
+ !!& info)
+
+ !write(*,*) ppmpoisson%fldzc1
+ #ifdef __NOPE
+ if (ppm_rank .EQ. trank) THEN
+ write(*,*)
+ write(*,*)
+ !ppmpoisson%fldzc1 = 0.0_MK
+ write(*,*) 'ifftx12', ppm_rank
+ DO isub=1,ppmpoisson%nsublistz
+ isubl=ppmpoisson%isublistz(isub)
+ DO k=1,ppmpoisson%ndataz(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataz(1,isubl)-1
+ DO j=1,ppmpoisson%ndataz(2,isubl)-1
+ write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldzc1(1,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ write(*,*) 'iffty12', ppm_rank
+ DO isub=1,ppmpoisson%nsublistz
+ isubl=ppmpoisson%isublistz(isub)
+ DO k=1,ppmpoisson%ndataz(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataz(1,isubl)-1
+ DO j=1,ppmpoisson%ndataz(2,isubl)-1
+ write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldzc1(2,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ write(*,*) 'ifftz12', ppm_rank
+ DO isub=1,ppmpoisson%nsublistz
+ isubl=ppmpoisson%isublistz(isub)
+ DO k=1,ppmpoisson%ndataz(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataz(1,isubl)-1
+ DO j=1,ppmpoisson%ndataz(2,isubl)-1
+ write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldzc1(3,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ ENDIF
+ #endif
+ !AOK
+
+ !!IF (isubl .EQ. 1) THEN
+ !!ppmpoisson%fldzc1 = CMPLX(-1.0_MK,-1.0_MK,MK)
+ !!ELSE
+ !!ppmpoisson%fldzc1 = CMPLX(8.0_MK,8.0_MK,MK)
+ !!DO isub=1,ppmpoisson%nsublistz
+ !!isubl=ppmpoisson%isublistz(isub)
+ !!DO k=1,ppmpoisson%ndataz(3,isubl)
+ !!DO j=1,ppmpoisson%ndataz(2,isubl)
+ !!DO i=1,ppmpoisson%ndataz(1,isubl)
+ !!ppmpoisson%fldzc1(1,i,j,k,isub) = REAL(k,MK)
+ !!END DO
+ !!END DO
+ !!END DO
+ !!END DO
+ !!ENDIF
+ !-----------------------------------------------------------------------
+ ! Map back to slabs (XY)
+ !-----------------------------------------------------------------------
+ !Initialise
+ CALL ppm_map_field_global(&
+ & ppmpoisson%topoidz, &
+ & ppmpoisson%topoidxy, &
+ & ppmpoisson%meshidz, &
+ & ppmpoisson%meshidxyc,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to initialise field mapping.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Push the data
+ CALL ppm_map_field_push(&
+ & ppmpoisson%topoidz, &
+ & ppmpoisson%meshidz,ppmpoisson%fldzc1,__NCOM,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to push vector field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Send
+ CALL ppm_map_field_send(info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to send field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !!ppmpoisson%fldxyc = CMPLX(-1.0_MK,-1.0_MK,MK)
+ !Retrieve
+ CALL ppm_map_field_pop(&
+ & ppmpoisson%topoidxy, &
+ & ppmpoisson%meshidxyc,ppmpoisson%fldxyc, &
+ & __NCOM,__ZEROSI,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !!write(*,*) 'ifftz slab'
+ !!DO isub=1,ppmpoisson%nsublistxy
+ !!isubl=ppmpoisson%isublistxy(isub)
+ !!IF (isubl.EQ.2) THEN
+ !!DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ !!write(*,*) 'z',k
+ !!DO i=1,ppmpoisson%ndataxy(1,isubl)-1
+ !!DO j=1,ppmpoisson%ndataxy(2,isubl)-1
+ !!write(*,'(A,E12.4,E12.4,A,$)') '(',ppmpoisson%fldxyc(1,i,j,k,isub),')'
+ !!END DO
+ !!write(*,*)
+ !!END DO
+ !!END DO
+ !!ENDIF
+ !!END DO
+ !-----------------------------------------------------------------------
+ ! IFFT (XY) use the non-reduced topology
+ !-----------------------------------------------------------------------
+ CALL ppm_fft_execute_2d(ppmpoisson%topoidxy,&
+ & ppmpoisson%meshidxy, ppmpoisson%planbxy, &
+ & ppmpoisson%fldxyc, ppmpoisson%fldxyr, &
+ & info)
+
+ !!CALL ppm_fft_normalize(ppmpoisson%topoidxy,&
+ !!& ppmpoisson%meshidxy, ppmpoisson%planbxy, __ZEROSI,&
+ !!& ppmpoisson%fldxyr, &
+ !!& info)
+
+ #ifdef __NOPE
+ if (ppm_rank .EQ. trank) THEN
+ write(*,*)
+ write(*,*)
+ write(*,*) 'ifftx', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)-1
+ write(*,'(E12.4,$)') ppmpoisson%fldxyr(1,i,j,k,isub)
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ write(*,*) 'iffty', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)-1
+ write(*,'(E12.4,$)') ppmpoisson%fldxyr(2,i,j,k,isub)
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ write(*,*) 'ifftz', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)-1
+ write(*,'(E12.4,$)') ppmpoisson%fldxyr(3,i,j,k,isub)
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+ ENDIF
+ #endif
+
+ !-----------------------------------------------------------------------
+ ! Map back to standard topology (XYZ)
+ !-----------------------------------------------------------------------
+ !Initialise
+ CALL ppm_map_field_global(&
+ & ppmpoisson%topoidxy, &
+ & topoid, &
+ & ppmpoisson%meshidxy, &
+ & meshid,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to initialise field mapping.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Push the data
+ CALL ppm_map_field_push(&
+ & ppmpoisson%topoidxy, &
+ & ppmpoisson%meshidxy,ppmpoisson%fldxyr,__NCOM,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to push vector field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Send
+ CALL ppm_map_field_send(info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to send field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! FINAL RETRIEVE - Here we do different things depending on the task
+ ! i.e. the receiver varies
+ !-------------------------------------------------------------------------
+ IF (ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_fd2 .AND.&
+ presentcase .EQ. ppm_poisson_grn_pois_per) THEN
+ CALL ppm_map_field_pop(&
+ & topoid, &
+ & meshid,ppmpoisson%drv_vr, &
+ & __NCOM,gstw,info) !!! gst
+ !-------------------------------------------------------------------------
+ ! Ghost the temporary array for derivatives (drv_vr)
+ !-------------------------------------------------------------------------
+ CALL ppm_map_field_ghost_get(topoid,meshid,gstw,info)
+ CALL ppm_map_field_push(topoid,meshid,ppmpoisson%drv_vr,__NCOM,info)
+ CALL ppm_map_field_send(info)
+ CALL ppm_map_field_pop(topoid,meshid,ppmpoisson%drv_vr,__NCOM,gstw,info)
+
+ ELSE
+ CALL ppm_map_field_pop(&
+ & topoid, &
+ & meshid,fieldout, &
+ & __NCOM,gstw,info) !!! gst
+ ENDIF
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
+ GOTO 9999
+ ENDIF
+
+
+ !-------------------------------------------------------------------------
+ ! Optionally do derivatives
+ ! Perhaps make ppm_poisson_fd take _none as argument. Then maybe no
+ ! if-statement is required
+ !-------------------------------------------------------------------------
+ IF (ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_fd2 .AND.&
+ presentcase .EQ. ppm_poisson_grn_pois_per) THEN
+ CALL ppm_poisson_fd(topoid,meshid,ppmpoisson%drv_vr,fieldout,&
+ & ppm_poisson_drv_curl_fd2,info)
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Finally ghost the velocity/stream function field before returning it
+ !-------------------------------------------------------------------------
+ CALL ppm_map_field_ghost_get(topoid,meshid,gstw,info)
+ CALL ppm_map_field_push(topoid,meshid,fieldout,__NCOM,info)
+ CALL ppm_map_field_send(info)
+ CALL ppm_map_field_pop(topoid,meshid,fieldout,__NCOM,gstw,info)
+
+
+ !-------------------------------------------------------------------------
+ ! Perhaps allocate (and deallocate) arrays !@
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_poisson_solve',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
+
diff --git a/src/ppm_define.h b/src/ppm_define.h
index e69de29..aa0a6dc 100644
--- a/src/ppm_define.h
+++ b/src/ppm_define.h
@@ -0,0 +1,2 @@
+#define __MPI
+#define __Linux
diff --git a/src/ppm_module_fft.f90 b/src/ppm_module_fft.f90
new file mode 100644
index 0000000..4ddca4e
--- /dev/null
+++ b/src/ppm_module_fft.f90
@@ -0,0 +1,154 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_module_fft
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! This modules contains routines for creating FFTW plans and executing them
+ ! Plan routines are required for
+ ! z means fft in x (1D)
+ ! xy means fft in xy (2D)
+ ! xyz means fft in xyz (3D)
+ ! sca means scalar
+ ! vec means vector
+ ! f means forward fft
+ ! b means backward fft
+ ! c means complex
+ ! r means real
+ ! s means single precision
+ ! d means double precision
+ !
+ ! e.g. 3d_vec_fr2c_xy_d: 3D vector array forward transform from
+ ! real to complex and ffts in x and y directions
+ ! (2d). For double data
+ !
+ ! The following variants have been implemented
+ ! 3d_vec_fc2c_z_d
+ ! 3d_vec_bc2c_z_d
+ ! 3d_vec_fr2c_xy_d
+ ! 3d_vec_bc2r_xy_d
+ !
+ ! a normalization routine exists for debugging purposes
+ !
+ ! sofar all fftw calls are to double precision routines!
+ !-------------------------------------------------------------------------
+#define __SINGLE 0
+#define __DOUBLE 1
+
+ MODULE ppm_module_fft
+ !-------------------------------------------------------------------------
+ ! 1D transforms
+ !-------------------------------------------------------------------------
+ INTERFACE ppm_fft_forward_1d
+ MODULE PROCEDURE ppm_fft_plan_3d_vec_fc2c_z_s
+ MODULE PROCEDURE ppm_fft_plan_3d_vec_fc2c_z_d
+ END INTERFACE
+
+ INTERFACE ppm_fft_backward_1d
+ MODULE PROCEDURE ppm_fft_plan_3d_vec_bc2c_z_s
+ MODULE PROCEDURE ppm_fft_plan_3d_vec_bc2c_z_d
+ END INTERFACE
+
+ INTERFACE ppm_fft_execute_1d
+ MODULE PROCEDURE ppm_fft_exec_3d_vec_c2c_z_s
+ MODULE PROCEDURE ppm_fft_exec_3d_vec_c2c_z_d
+ END INTERFACE
+
+ !-------------------------------------------------------------------------
+ ! 2D transforms
+ !-------------------------------------------------------------------------
+ INTERFACE ppm_fft_forward_2d
+ MODULE PROCEDURE ppm_fft_plan_3d_vec_fr2c_xy_s
+ MODULE PROCEDURE ppm_fft_plan_3d_vec_fr2c_xy_d
+ END INTERFACE
+
+ INTERFACE ppm_fft_backward_2d
+ MODULE PROCEDURE ppm_fft_plan_3d_vec_bc2r_xy_s
+ MODULE PROCEDURE ppm_fft_plan_3d_vec_bc2r_xy_d
+ END INTERFACE
+
+ INTERFACE ppm_fft_execute_2d
+ MODULE PROCEDURE ppm_fft_exec_3d_vec_fr2c_xy_s
+ MODULE PROCEDURE ppm_fft_exec_3d_vec_bc2r_xy_s
+ MODULE PROCEDURE ppm_fft_exec_3d_vec_fr2c_xy_d
+ MODULE PROCEDURE ppm_fft_exec_3d_vec_bc2r_xy_d
+ END INTERFACE
+
+ !-------------------------------------------------------------------------
+ ! Normalization
+ !-------------------------------------------------------------------------
+ INTERFACE ppm_fft_normalize
+ MODULE PROCEDURE ppm_fft_normalize_rs
+ MODULE PROCEDURE ppm_fft_normalize_cs
+ MODULE PROCEDURE ppm_fft_normalize_rd
+ MODULE PROCEDURE ppm_fft_normalize_cd
+ END INTERFACE
+
+ !-------------------------------------------------------------------------
+ ! PPM FFT plan type
+ !-------------------------------------------------------------------------
+ !!! Type containing the FFTW plan and its settings
+ TYPE ppm_fft_plan
+ !!!array of plan pointers, index for subs
+ INTEGER*8,DIMENSION(:),POINTER :: plan => NULL()
+ !!!the dimension of the FFT (1D/2D/3D)
+ INTEGER :: rank
+ !!!number of points along each direction of the piece to be transformed
+ !!!index is for rank and subs
+ INTEGER,DIMENSION(:,:),POINTER :: nx
+ !!!the direction of the transform (forward/backward)
+ INTEGER :: sign
+ !!!the method to setup the optimal plan
+ INTEGER :: flag
+ !!!the number of components to transform
+ INTEGER :: howmany
+ !!!the size of the input array, index is for rank
+ INTEGER,DIMENSION(:),POINTER :: inembed
+ !!!the size of the output array, index is for rank
+ INTEGER,DIMENSION(:),POINTER :: onembed
+ !!!istride tells how same component data points are spaced in memory
+ INTEGER :: istride
+ INTEGER :: ostride
+ !!!idist tells how multiple arrays are spaced. I.e. a memory offset
+ INTEGER :: idist
+ INTEGER :: odist
+ END TYPE ppm_fft_plan
+
+ CONTAINS
+#define __KIND __SINGLE
+
+ !FORWARD TRANSFORMS - PLAN
+#include "fft/ppm_fft_plan_3d_vec_fc2c_z.f90"
+#include "fft/ppm_fft_plan_3d_vec_fr2c_xy.f90"
+ !BACKWARD TRANSFORMS - PLAN
+#include "fft/ppm_fft_plan_3d_vec_bc2c_z.f90"
+#include "fft/ppm_fft_plan_3d_vec_bc2r_xy.f90"
+ !EXECUTION OF PLANS
+#include "fft/ppm_fft_exec_3d_vec_c2c_z.f90"
+#include "fft/ppm_fft_exec_3d_vec_fr2c_xy.f90"
+#include "fft/ppm_fft_exec_3d_vec_bc2r_xy.f90"
+ !NORMALIZATION
+#include "fft/ppm_fft_normalize_r.f90"
+#include "fft/ppm_fft_normalize_c.f90"
+
+#undef __KIND
+#define __KIND __DOUBLE
+
+ !FORWARD TRANSFORMS - PLAN
+#include "fft/ppm_fft_plan_3d_vec_fc2c_z.f90"
+#include "fft/ppm_fft_plan_3d_vec_fr2c_xy.f90"
+ !BACKWARD TRANSFORMS - PLAN
+#include "fft/ppm_fft_plan_3d_vec_bc2c_z.f90"
+#include "fft/ppm_fft_plan_3d_vec_bc2r_xy.f90"
+ !EXECUTION OF PLANS
+#include "fft/ppm_fft_exec_3d_vec_c2c_z.f90"
+#include "fft/ppm_fft_exec_3d_vec_fr2c_xy.f90"
+#include "fft/ppm_fft_exec_3d_vec_bc2r_xy.f90"
+ !NORMALIZATION
+#include "fft/ppm_fft_normalize_r.f90"
+#include "fft/ppm_fft_normalize_c.f90"
+
+#undef __KIND
+ END MODULE ppm_module_fft
+
+
diff --git a/src/ppm_module_poisson.f90 b/src/ppm_module_poisson.f90
new file mode 100644
index 0000000..2c7d8e8
--- /dev/null
+++ b/src/ppm_module_poisson.f90
@@ -0,0 +1,195 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_module_poisson
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
+ ! Notes on the init routine:
+ ! Arguments:
+ ! topoid
+ ! meshid
+ ! ppmpoisson
+ ! fieldin
+ ! fieldout
+ ! greens function: custom array, custom function, integer to build-in
+ ! info - return variable
+ ! periodic/freespace flag
+ ! real/fourier greens function (optional, only for custom array)
+ ! derivatives (optional): none,curl,gradient,laplace, - spectral/FD
+ ! possibly the option to de- and reallocate arrays of slabs/pencils
+ !
+ !
+ ! Add finalize routine
+ !-------------------------------------------------------------------------
+#define __SINGLE 0
+#define __DOUBLE 1
+
+#define __DIM 3
+
+ MODULE ppm_module_poisson
+
+ USE ppm_module_fft
+ USE ppm_module_substart
+ USE ppm_module_substop
+ USE ppm_module_write
+ USE ppm_module_data,ONLY:ppm_rank,ppm_kind_single,ppm_kind_double,&
+ &ppm_t_topo,ppm_t_equi_mesh,&
+ &ppm_param_assign_internal,ppm_param_bcdef_periodic,&
+ &ppm_param_bcdef_freespace,&
+ &ppm_param_decomp_xy_slab,ppm_param_decomp_zpencil
+
+ !-------------------------------------------------------------------------
+ ! PPM Poisson type
+ !-------------------------------------------------------------------------
+ !!! Type containing the FFTW plan and its settings
+ TYPE ppm_poisson_plan
+ INTEGER :: derivatives
+ INTEGER :: case
+
+
+ INTEGER :: topoidxy
+ !!!Topology id of xy topology
+ INTEGER :: meshidxy
+ !!!mesh id for xy mesh
+ REAL(ppm_kind_double), DIMENSION(:), POINTER :: costxy=>NULL()
+ !!!sub cost for xy topology
+ INTEGER , DIMENSION(:,:),POINTER :: istartxy=>NULL()
+ !!!sub index start for xy topology
+ INTEGER, DIMENSION(__DIM) :: nmxy
+ !!!global number of grid points for xy topology
+ INTEGER , DIMENSION(:,:),POINTER :: ndataxy=>NULL()
+ !!!sub no. grid cells for xy topology
+ INTEGER,DIMENSION(:),POINTER :: isublistxy
+ !!!list of sub domains on current CPU on the xy topology
+ INTEGER :: nsublistxy
+ !!!number of sub domains for current CPU on the xy topology
+ REAL(ppm_kind_double),DIMENSION(:,:,:,:,:),POINTER :: fldxyr=>NULL()
+ !!!real slab field (xy topology)
+ COMPLEX(ppm_kind_double),DIMENSION(:,:,:,:,:),POINTER:: fldxyc=>NULL()
+ !!!complex slab field (xy topology)
+
+ INTEGER :: topoidxyc
+ !!!Topology id of complex xy topology
+ INTEGER :: meshidxyc
+ !!!mesh id for comple xy mesh
+ REAL(ppm_kind_double), DIMENSION(:), POINTER :: costxyc=>NULL()
+ !!!sub cost for xy topology
+ INTEGER , DIMENSION(:,:),POINTER :: istartxyc=>NULL()
+ !!!sub index start for complex xy topology
+ INTEGER, DIMENSION(__DIM) :: nmxyc
+ !!!global number of grid points for complex xy topology
+ INTEGER , DIMENSION(:,:),POINTER :: ndataxyc=>NULL()
+ !!!sub no. grid cells for complex xy topology
+ INTEGER,DIMENSION(:),POINTER :: isublistxyc
+ !!!list of sub domains on current CPU on the complex xy topology
+ INTEGER :: nsublistxyc
+ !!!number of sub domains for current CPU on the complex xy topology
+ TYPE(ppm_fft_plan) :: planfxy
+ !!!fft plans for r2c complex xy slab
+ TYPE(ppm_fft_plan) :: planbxy
+ !!!fft plans for c2r complex xy slab
+
+ TYPE(ppm_fft_plan) :: planfz
+ !!!fft plans for forward c2c z pencil
+ TYPE(ppm_fft_plan) :: planbz
+ !!!fft plans for backward c2c z pencil
+ INTEGER :: topoidz
+ !!!Topology id of z topology
+ INTEGER :: meshidz
+ !!!mesh id for z mesh
+ REAL(ppm_kind_double), DIMENSION(:), POINTER :: costz =>NULL()
+ !!!sub cost for z topology
+ INTEGER , DIMENSION(:,:),POINTER :: istartz=>NULL()
+ !!!sub index start for z topology
+ INTEGER, DIMENSION(__DIM) :: nmz
+ !!!global number of grid points for z topology
+ INTEGER , DIMENSION(:,:),POINTER :: ndataz=>NULL()
+ !!!sub no. grid cells for z topology
+ INTEGER,DIMENSION(:),POINTER :: isublistz
+ !!!list of sub domains on current CPU on the z topology
+ INTEGER :: nsublistz
+ !!!number of sub domains for current CPU on the z topology
+ COMPLEX(ppm_kind_double),DIMENSION(:,:,:,:,:),POINTER:: fldzc1=>NULL()
+ !!!complex pencil field 1 (z topology)
+ COMPLEX(ppm_kind_double),DIMENSION(:,:,:,:,:),POINTER:: fldzc2=>NULL()
+ !!!complex pencil field 2 (z topology)
+
+ INTEGER :: topoidfr
+ !!!Topology id of freespace topology
+ INTEGER :: meshidfr
+ !!!mesh id for freespace mesh
+ REAL(ppm_kind_double), DIMENSION(:), POINTER :: costfr=>NULL()
+ !!!sub cost
+ INTEGER , DIMENSION(:,:),POINTER :: istartfr=>NULL()
+ !!!sub index start
+ INTEGER, DIMENSION(__DIM) :: nmfr
+ !!!global number of grid points
+ INTEGER , DIMENSION(:,:),POINTER :: ndatafr=>NULL()
+ !!!sub no. grid cells
+ INTEGER,DIMENSION(:),POINTER :: isublistfr
+ !!!list of sub domains on current CPU on the freespace topology
+ INTEGER :: nsublistfr
+ !!!number of sub domains for current CPU on the freespace topology
+
+ !These field are only allocated as necessary:
+ REAL(ppm_kind_double),DIMENSION(:,:,:,:),POINTER :: fldfrs=>NULL()
+ !!!dummy array for the right hand side, for free space, scalar fields
+ REAL(ppm_kind_double),DIMENSION(:,:,:,:,:),POINTER :: fldfrv=>NULL()
+ !!!dummy array for the right hand side, for free space, vector fields
+ REAL(ppm_kind_double),DIMENSION(:,:,:,:),POINTER :: fldgrnr=>NULL()
+ !!!real Greens field, z-pencils, scalar
+ COMPLEX(ppm_kind_double),DIMENSION(:,:,:,:,:),POINTER :: fldgrnc=>NULL() !@tmp
+ !!!complex Greens field, z-pencils, scalar
+ REAL(ppm_kind_double),DIMENSION(:,:,:,:,:),POINTER :: drv_vr=>NULL()
+ !!!dummy array for the right hand side, for free space, vector fields
+
+ ! COMPLEX(ppm_kind_double),DIMENSION(:,:,:,:,:),POINTER:: fldgrnc3=>NULL()
+ ! !!!complex Greens 3 component vector field intended for spectral
+ ! !!!derivatives, z-pencils, scalar
+ END TYPE ppm_poisson_plan
+
+ INTEGER,PARAMETER :: ppm_poisson_drv_none =0
+ INTEGER,PARAMETER :: ppm_poisson_drv_curl_sp =1
+ INTEGER,PARAMETER :: ppm_poisson_drv_grad_sp =2
+ INTEGER,PARAMETER :: ppm_poisson_drv_lapl_sp =3
+ INTEGER,PARAMETER :: ppm_poisson_drv_div_sp =4
+ INTEGER,PARAMETER :: ppm_poisson_drv_curl_fd2=11
+ INTEGER,PARAMETER :: ppm_poisson_drv_grad_fd2=12
+ INTEGER,PARAMETER :: ppm_poisson_drv_lapl_fd2=13
+ INTEGER,PARAMETER :: ppm_poisson_drv_div_fd2 =14
+ INTEGER,PARAMETER :: ppm_poisson_drv_curl_fd4=21
+ INTEGER,PARAMETER :: ppm_poisson_drv_grad_fd4=22
+ INTEGER,PARAMETER :: ppm_poisson_drv_lapl_fd4=23
+ INTEGER,PARAMETER :: ppm_poisson_drv_div_fd4 =24
+
+ INTEGER,PARAMETER :: ppm_poisson_grn_pois_per =1
+ INTEGER,PARAMETER :: ppm_poisson_grn_pois_fre =2
+ INTEGER,PARAMETER :: ppm_poisson_grn_reprojec =3
+
+ INTERFACE ppm_poisson_init
+ MODULE PROCEDURE ppm_poisson_init_predef
+ END INTERFACE
+
+ INTERFACE ppm_poisson_solve
+ MODULE PROCEDURE ppm_poisson_solve
+ END INTERFACE
+
+ INTERFACE ppm_poisson_fd
+ MODULE PROCEDURE ppm_poisson_fd
+ END INTERFACE
+
+ CONTAINS
+#define __KIND __SINGLE
+
+#include "poisson/ppm_poisson_init_predef.f90"
+#include "poisson/ppm_poisson_solve.f90"
+#include "poisson/ppm_poisson_fd.f90"
+
+#undef __KIND
+#define __KIND __DOUBLE
+
+
+#undef __KIND
+ END MODULE ppm_module_poisson
+
+
--
GitLab