From 84c9c96c49be4d78ad34dda5feed83467ae983e6 Mon Sep 17 00:00:00 2001
From: jtra <jtra@7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9>
Date: Wed, 1 Jun 2011 07:49:45 +0000
Subject: [PATCH] working
git-svn-id: https://ppm.inf.ethz.ch/svn/ppmnumerics/branches/ngtopo/libppmnumerics@884 7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9
---
src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f | 119 +++
src/fft/ppm_fft_exec_3d_vec_c2c_z.f | 124 +++
src/fft/ppm_fft_exec_3d_vec_fr2c_xy.f | 103 +++
src/fft/ppm_fft_normalize_c.f | 127 ++++
src/fft/ppm_fft_normalize_r.f | 128 ++++
src/fft/ppm_fft_plan_3d_vec_bc2c_z.f | 149 ++++
src/fft/ppm_fft_plan_3d_vec_bc2r_xy.f | 154 ++++
src/fft/ppm_fft_plan_3d_vec_fc2c_z.f | 149 ++++
src/fft/ppm_fft_plan_3d_vec_fr2c_xy.f | 153 ++++
src/poisson/ppm_poisson_extrapolateghost.f | 257 +++++++
src/poisson/ppm_poisson_fd.f | 147 ++++
src/poisson/ppm_poisson_init_predef.f | 637 ++++++++++++++++
src/poisson/ppm_poisson_solve.f | 837 +++++++++++++++++++++
13 files changed, 3084 insertions(+)
create mode 100644 src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f
create mode 100644 src/fft/ppm_fft_exec_3d_vec_c2c_z.f
create mode 100644 src/fft/ppm_fft_exec_3d_vec_fr2c_xy.f
create mode 100644 src/fft/ppm_fft_normalize_c.f
create mode 100644 src/fft/ppm_fft_normalize_r.f
create mode 100644 src/fft/ppm_fft_plan_3d_vec_bc2c_z.f
create mode 100644 src/fft/ppm_fft_plan_3d_vec_bc2r_xy.f
create mode 100644 src/fft/ppm_fft_plan_3d_vec_fc2c_z.f
create mode 100644 src/fft/ppm_fft_plan_3d_vec_fr2c_xy.f
create mode 100644 src/poisson/ppm_poisson_extrapolateghost.f
create mode 100644 src/poisson/ppm_poisson_fd.f
create mode 100644 src/poisson/ppm_poisson_init_predef.f
create mode 100644 src/poisson/ppm_poisson_solve.f
diff --git a/src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f b/src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f
new file mode 100644
index 0000000..de26a49
--- /dev/null
+++ b/src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f
@@ -0,0 +1,119 @@
+ !-------------------------------------------------------------------------
+ ! 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))
+ IF (topology%bcdef(1) .EQ. ppm_param_bcdef_periodic) THEN !vertex
+ !-------------------------------------------------------------------
+ ! 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 IF
+ 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.f b/src/fft/ppm_fft_exec_3d_vec_c2c_z.f
new file mode 100644
index 0000000..ca01f76
--- /dev/null
+++ b/src/fft/ppm_fft_exec_3d_vec_c2c_z.f
@@ -0,0 +1,124 @@
+ !-------------------------------------------------------------------------
+ ! 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 (topology%bcdef(3) .EQ. ppm_param_bcdef_periodic) THEN !vertex
+ 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
+ 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.f b/src/fft/ppm_fft_exec_3d_vec_fr2c_xy.f
new file mode 100644
index 0000000..82ea0e3
--- /dev/null
+++ b/src/fft/ppm_fft_exec_3d_vec_fr2c_xy.f
@@ -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.f b/src/fft/ppm_fft_normalize_c.f
new file mode 100644
index 0000000..d8f228f
--- /dev/null
+++ b/src/fft/ppm_fft_normalize_c.f
@@ -0,0 +1,127 @@
+ !-------------------------------------------------------------------------
+ ! 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
+ IF (topology%bcdef(3) .EQ. ppm_param_bcdef_periodic) THEN !vertex
+ fac = 1.0_MK/REAL((mesh%nm(3)-1),MK)
+ ELSE
+ fac = 1.0_MK/REAL((mesh%nm(3)),MK)
+ ENDIF
+ ELSE IF (ppmplan%rank .EQ. 2) THEN
+ IF (topology%bcdef(1) .EQ. ppm_param_bcdef_periodic) THEN !vertex
+ fac = 1.0_MK/REAL((mesh%nm(1)-1)*(mesh%nm(2)-1),MK)
+ ELSE
+ fac = 1.0_MK/REAL((mesh%nm(1) )*(mesh%nm(2) ),MK)
+ ENDIF
+ 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.f b/src/fft/ppm_fft_normalize_r.f
new file mode 100644
index 0000000..425ec12
--- /dev/null
+++ b/src/fft/ppm_fft_normalize_r.f
@@ -0,0 +1,128 @@
+ !-------------------------------------------------------------------------
+ ! 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
+ IF (topology%bcdef(3) .EQ. ppm_param_bcdef_periodic) THEN !vertex
+ fac = 1.0_MK/REAL((mesh%nm(3)-1),MK)
+ ELSE
+ fac = 1.0_MK/REAL((mesh%nm(3)),MK)
+ ENDIF
+ ELSE IF (ppmplan%rank .EQ. 2) THEN
+ IF (topology%bcdef(1) .EQ. ppm_param_bcdef_periodic) THEN !vertex
+ fac = 1.0_MK/REAL((mesh%nm(1)-1)*(mesh%nm(2)-1),MK)
+ ELSE
+ fac = 1.0_MK/REAL((mesh%nm(1) )*(mesh%nm(2) ),MK)
+ ENDIF
+ 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.f b/src/fft/ppm_fft_plan_3d_vec_bc2c_z.f
new file mode 100644
index 0000000..9dd2ba2
--- /dev/null
+++ b/src/fft/ppm_fft_plan_3d_vec_bc2c_z.f
@@ -0,0 +1,149 @@
+ !-------------------------------------------------------------------------
+ ! 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
+ IF (topology%bcdef(3) .EQ. ppm_param_bcdef_periodic) THEN !vertex
+ ppmplan%nx(1,isub) = mesh%nnodes(3,isubl)-1
+ ELSE
+ ppmplan%nx(1,isub) = mesh%nnodes(3,isubl)
+ ENDIF
+
+ 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.f b/src/fft/ppm_fft_plan_3d_vec_bc2r_xy.f
new file mode 100644
index 0000000..f83a32e
--- /dev/null
+++ b/src/fft/ppm_fft_plan_3d_vec_bc2r_xy.f
@@ -0,0 +1,154 @@
+ !-------------------------------------------------------------------------
+ ! 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
+ IF (topology%bcdef(1) .EQ. ppm_param_bcdef_periodic) THEN !vertex
+ ppmplan%nx(1,isub) = mesh%nnodes(1,isubl)-1
+ ppmplan%nx(2,isub) = mesh%nnodes(2,isubl)-1
+ ELSE
+ ppmplan%nx(1,isub) = mesh%nnodes(1,isubl)
+ ppmplan%nx(2,isub) = mesh%nnodes(2,isubl)
+ ENDIF
+
+ 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.f b/src/fft/ppm_fft_plan_3d_vec_fc2c_z.f
new file mode 100644
index 0000000..2625176
--- /dev/null
+++ b/src/fft/ppm_fft_plan_3d_vec_fc2c_z.f
@@ -0,0 +1,149 @@
+ !-------------------------------------------------------------------------
+ ! 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
+ IF (topology%bcdef(3) .EQ. ppm_param_bcdef_periodic) THEN !vertex
+ ppmplan%nx(1,isub) = mesh%nnodes(3,isubl)-1
+ ELSE
+ ppmplan%nx(1,isub) = mesh%nnodes(3,isubl)
+ ENDIF
+
+ 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.f b/src/fft/ppm_fft_plan_3d_vec_fr2c_xy.f
new file mode 100644
index 0000000..d1eab5b
--- /dev/null
+++ b/src/fft/ppm_fft_plan_3d_vec_fr2c_xy.f
@@ -0,0 +1,153 @@
+ !-------------------------------------------------------------------------
+ ! 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
+ IF (topology%bcdef(1) .EQ. ppm_param_bcdef_periodic) THEN !vertex
+ ppmplan%nx(1,isub) = mesh%nnodes(1,isubl)-1
+ ppmplan%nx(2,isub) = mesh%nnodes(2,isubl)-1
+ ELSE
+ ppmplan%nx(1,isub) = mesh%nnodes(1,isubl)
+ ppmplan%nx(2,isub) = mesh%nnodes(2,isubl)
+ ENDIF
+
+ 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_extrapolateghost.f b/src/poisson/ppm_poisson_extrapolateghost.f
new file mode 100644
index 0000000..c2cc564
--- /dev/null
+++ b/src/poisson/ppm_poisson_extrapolateghost.f
@@ -0,0 +1,257 @@
+ !-------------------------------------------------------------------------
+ ! ppm_poisson_extrapolateghost.f90
+ !-------------------------------------------------------------------------
+!!#define __ROUTINE ppm_poisson_extrapolateghost_vr
+!!#define __DIM 3
+!!#define __NCOM 3
+!!#define __ZEROSI (/0,0,0/)
+ SUBROUTINE __ROUTINE(topoid,meshid,field,nextra,nbase,gstw,info)
+
+ USE ppm_module_topo_get
+
+
+ IMPLICIT NONE
+ !-------------------------------------------------------------------------
+ ! Arguments
+ !-------------------------------------------------------------------------
+ INTEGER, INTENT(IN) :: topoid
+ !!! Topology id of the field
+ INTEGER, INTENT(IN) :: meshid
+ !!! Mesh id of the field
+ REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: field
+ !!! Field to extrapolate ghost layer
+ INTEGER, INTENT(IN) :: nextra
+ !!! Number of points to extrapolate into the ghostlayer
+ INTEGER, INTENT(IN) :: nbase
+ !!! Number of points to base the extrapolation on
+ INTEGER,DIMENSION(__DIM),INTENT(IN) :: gstw
+ !!! Width of the ghotslayer
+ INTEGER, INTENT(OUT) :: info
+ !!! Return state
+
+ !-------------------------------------------------------------------------
+ ! Local variables
+ !-------------------------------------------------------------------------
+ INTEGER,PARAMETER :: MK = __PREC
+ REAL(__PREC) :: t0
+ TYPE(ppm_t_topo),POINTER :: topology
+ TYPE(ppm_t_equi_mesh) :: mesh
+ INTEGER :: isub,isubl
+ INTEGER :: i,j,k,iextra,ibase
+ REAL(__PREC),DIMENSION(:,:),POINTER :: coeff
+ REAL(__PREC),DIMENSION(__DIM) :: tmpbuf
+
+
+ !-------------------------------------------------------------------------
+ ! Initialise routine
+ !-------------------------------------------------------------------------
+ CALL substart('ppm_poisson_extrapolateghost',t0,info)
+
+
+ !-------------------------------------------------------------------------
+ ! Compare the number of points to extrapolate to the ghost layer width
+ !-------------------------------------------------------------------------
+ IF (iextra .GT. gstw(1) .OR. &
+ & iextra .GT. gstw(2) .OR. &
+ & iextra .GT. gstw(3)) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_extrapolateghost',&
+ & 'The points to extrapolate exceeds the ghost layer.',info)
+ info = -1
+ GOTO 9999
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ !@ Check what has been implemented in this routine so far
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Determine weights
+ !-------------------------------------------------------------------------
+ ALLOCATE(coeff(nbase,nextra))
+ IF (nbase .EQ. 4) THEN
+ IF (nextra .GE. 1) THEN
+ coeff(:,1) = (/4.0_MK,-6.0_MK,4.0_MK,-1.0_MK/)
+ ENDIF
+ IF (nextra .GE. 2) THEN
+ coeff(:,2) = (/10.0_MK,-20.0_MK,15.0_MK,-4.0_MK/)
+ ENDIF
+ IF (nextra .GE. 3) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_extrapolateghost',&
+ & 'Extrapolation to more than two points has not been implemented.',info)
+ info = -1
+ GOTO 9999
+ ENDIF
+ ELSE
+ CALL ppm_write(ppm_rank,'ppm_poisson_extrapolateghost',&
+ & 'Only extrapolation based on 4 points has been implemented.',info)
+ info = -1
+ GOTO 9999
+ ENDIF
+
+
+ !-------------------------------------------------------------------------
+ ! Get topology
+ !-------------------------------------------------------------------------
+ CALL ppm_topo_get(topoid,topology,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_extrapolateghost',&
+ & 'Failed to get topology.',isub)
+ GOTO 9999
+ ENDIF
+ mesh = topology%mesh(meshid)
+
+ !-------------------------------------------------------------------------
+ ! Extrapolate field into ghost layer
+ ! The indicies of subs_bc represent:
+ ! west,east(x),south,north(y),bottom,top(z)
+ !@ Some more unrolling here would be nice
+ !-------------------------------------------------------------------------
+ DO isub=1,topology%nsublist
+ isubl=topology%isublist(isub)
+ !West (-x)
+ IF (topology%subs_bc(1,isubl) .EQ. 1) THEN
+ DO k=1-gstw(3),mesh%nnodes(3,isubl)+gstw(3)
+ DO j=1-gstw(2),mesh%nnodes(2,isubl)+gstw(2)
+ !DO i=1-gstw(1),mesh%nnodes(1,isubl)+gstw(1)
+ i = 1
+ DO iextra=1,nextra
+ tmpbuf = 0.0_MK
+ DO ibase=0,nbase-1
+ tmpbuf(1) = tmpbuf(1) + &
+ & coeff(ibase+1,iextra)*field(1,i+ibase,j,k,isub)
+ tmpbuf(2) = tmpbuf(2) + &
+ & coeff(ibase+1,iextra)*field(2,i+ibase,j,k,isub)
+ tmpbuf(3) = tmpbuf(3) + &
+ & coeff(ibase+1,iextra)*field(3,i+ibase,j,k,isub)
+ END DO !ibase
+ field(1,i-iextra,j,k,isub) = tmpbuf(1)
+ field(2,i-iextra,j,k,isub) = tmpbuf(2)
+ field(3,i-iextra,j,k,isub) = tmpbuf(3)
+ END DO !iextra
+ ENDDO !j
+ ENDDO !k
+ ENDIF
+ !East (+x)
+ IF (topology%subs_bc(2,isubl) .EQ. 1) THEN
+ DO k=1-gstw(3),mesh%nnodes(3,isubl)+gstw(3)
+ DO j=1-gstw(2),mesh%nnodes(2,isubl)+gstw(2)
+ !DO i=1-gstw(1),mesh%nnodes(1,isubl)+gstw(1)
+ i = mesh%nnodes(1,isubl)
+ DO iextra=1,nextra
+ tmpbuf = 0.0_MK
+ DO ibase=0,nbase-1
+ tmpbuf(1) = tmpbuf(1) + &
+ & coeff(ibase+1,iextra)*field(1,i-ibase,j,k,isub)
+ tmpbuf(2) = tmpbuf(2) + &
+ & coeff(ibase+1,iextra)*field(2,i-ibase,j,k,isub)
+ tmpbuf(3) = tmpbuf(3) + &
+ & coeff(ibase+1,iextra)*field(3,i-ibase,j,k,isub)
+ END DO !ibase
+ field(1,i+iextra,j,k,isub) = tmpbuf(1)
+ field(2,i+iextra,j,k,isub) = tmpbuf(2)
+ field(3,i+iextra,j,k,isub) = tmpbuf(3)
+ END DO !iextra
+ ENDDO !j
+ ENDDO !k
+ ENDIF
+ !South (-y)
+ IF (topology%subs_bc(3,isubl) .EQ. 1) THEN
+ DO k=1-gstw(3),mesh%nnodes(3,isubl)+gstw(3)
+ !DO j=1-gstw(2),mesh%nnodes(2,isubl)+gstw(2)
+ DO i=1-gstw(1),mesh%nnodes(1,isubl)+gstw(1)
+ j = 1
+ DO iextra=1,nextra
+ tmpbuf = 0.0_MK
+ DO ibase=0,nbase-1
+ tmpbuf(1) = tmpbuf(1) + &
+ & coeff(ibase+1,iextra)*field(1,i,j+ibase,k,isub)
+ tmpbuf(2) = tmpbuf(2) + &
+ & coeff(ibase+1,iextra)*field(2,i,j+ibase,k,isub)
+ tmpbuf(3) = tmpbuf(3) + &
+ & coeff(ibase+1,iextra)*field(3,i,j+ibase,k,isub)
+ END DO !ibase
+ field(1,i,j-iextra,k,isub) = tmpbuf(1)
+ field(2,i,j-iextra,k,isub) = tmpbuf(2)
+ field(3,i,j-iextra,k,isub) = tmpbuf(3)
+ END DO !iextra
+ ENDDO !i
+ ENDDO !k
+ ENDIF
+ !North (+y)
+ IF (topology%subs_bc(4,isubl) .EQ. 1) THEN
+ DO k=1-gstw(3),mesh%nnodes(3,isubl)+gstw(3)
+ !DO j=1-gstw(2),mesh%nnodes(2,isubl)+gstw(2)
+ DO i=1-gstw(1),mesh%nnodes(1,isubl)+gstw(1)
+ j = mesh%nnodes(2,isubl)
+ DO iextra=1,nextra
+ tmpbuf = 0.0_MK
+ DO ibase=0,nbase-1
+ tmpbuf(1) = tmpbuf(1) + &
+ & coeff(ibase+1,iextra)*field(1,i,j-ibase,k,isub)
+ tmpbuf(2) = tmpbuf(2) + &
+ & coeff(ibase+1,iextra)*field(2,i,j-ibase,k,isub)
+ tmpbuf(3) = tmpbuf(3) + &
+ & coeff(ibase+1,iextra)*field(3,i,j-ibase,k,isub)
+ END DO !ibase
+ field(1,i,j+iextra,k,isub) = tmpbuf(1)
+ field(2,i,j+iextra,k,isub) = tmpbuf(2)
+ field(3,i,j+iextra,k,isub) = tmpbuf(3)
+ END DO !iextra
+ ENDDO !i
+ ENDDO !k
+ ENDIF
+ !Bottom (-z)
+ IF (topology%subs_bc(5,isubl) .EQ. 1) THEN
+ !DO k=1-gstw(3),mesh%nnodes(3,isubl)+gstw(3)
+ DO j=1-gstw(2),mesh%nnodes(2,isubl)+gstw(2)
+ DO i=1-gstw(1),mesh%nnodes(1,isubl)+gstw(1)
+ k = 1
+ DO iextra=1,nextra
+ tmpbuf = 0.0_MK
+ DO ibase=0,nbase-1
+ tmpbuf(1) = tmpbuf(1) + &
+ & coeff(ibase+1,iextra)*field(1,i,j,k+ibase,isub)
+ tmpbuf(2) = tmpbuf(2) + &
+ & coeff(ibase+1,iextra)*field(2,i,j,k+ibase,isub)
+ tmpbuf(3) = tmpbuf(3) + &
+ & coeff(ibase+1,iextra)*field(3,i,j,k+ibase,isub)
+ END DO !ibase
+ field(1,i,j,k-iextra,isub) = tmpbuf(1)
+ field(2,i,j,k-iextra,isub) = tmpbuf(2)
+ field(3,i,j,k-iextra,isub) = tmpbuf(3)
+ END DO !iextra
+ ENDDO !i
+ ENDDO !j
+ ENDIF
+ !Top (+z)
+ IF (topology%subs_bc(6,isubl) .EQ. 1) THEN
+ !DO k=1-gstw(3),mesh%nnodes(3,isubl)+gstw(3)
+ DO j=1-gstw(2),mesh%nnodes(2,isubl)+gstw(2)
+ DO i=1-gstw(1),mesh%nnodes(1,isubl)+gstw(1)
+ k = mesh%nnodes(3,isubl)
+ DO iextra=1,nextra
+ tmpbuf = 0.0_MK
+ DO ibase=0,nbase-1
+ tmpbuf(1) = tmpbuf(1) + &
+ & coeff(ibase+1,iextra)*field(1,i,j,k-ibase,isub)
+ tmpbuf(2) = tmpbuf(2) + &
+ & coeff(ibase+1,iextra)*field(2,i,j,k-ibase,isub)
+ tmpbuf(3) = tmpbuf(3) + &
+ & coeff(ibase+1,iextra)*field(3,i,j,k-ibase,isub)
+ END DO !ibase
+ field(1,i,j,k+iextra,isub) = tmpbuf(1)
+ field(2,i,j,k+iextra,isub) = tmpbuf(2)
+ field(3,i,j,k+iextra,isub) = tmpbuf(3)
+ END DO !iextra
+ ENDDO !i
+ ENDDO !j
+ ENDIF
+ ENDDO !isub
+
+
+ 9999 CONTINUE
+ CALL substop('ppm_poisson_extrapolateghost',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
+
diff --git a/src/poisson/ppm_poisson_fd.f b/src/poisson/ppm_poisson_fd.f
new file mode 100644
index 0000000..ebc7415
--- /dev/null
+++ b/src/poisson/ppm_poisson_fd.f
@@ -0,0 +1,147 @@
+ !-------------------------------------------------------------------------
+ ! ppm_poisson_fd.f90
+ !-------------------------------------------------------------------------
+ !@ TODO: Somewhere check if fieldin is equal to fieldout and give a warning
+ !-------------------------------------------------------------------------
+ !!#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)
+
+ !@ perhaps ensure that fieldin .NE. fieldout
+ !-------------------------------------------------------------------------
+ ! 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) !vertex
+ 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) = &
+ & facy*(fieldin(3,i ,j+1,k ,isub)- &
+ & fieldin(3,i ,j-1,k ,isub)) &
+ & -facz*(fieldin(2,i ,j ,k+1,isub)- &
+ & fieldin(2,i ,j ,k-1,isub))
+ fieldout(2,i,j,k,isub) = &
+ & facz*(fieldin(1,i ,j ,k+1,isub)- &
+ & fieldin(1,i ,j ,k-1,isub)) &
+ & -facx*(fieldin(3,i+1,j ,k ,isub)- &
+ & fieldin(3,i-1,j ,k ,isub))
+ fieldout(3,i,j,k,isub) = &
+ & facx*(fieldin(2,i+1,j ,k ,isub)- &
+ & fieldin(2,i-1,j ,k ,isub)) &
+ & -facy*(fieldin(1,i ,j+1,k ,isub)- &
+ & fieldin(1,i ,j-1,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) = &
+ & 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)) &
+ & -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))
+ fieldout(2,i,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)) &
+ & -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))
+ fieldout(3,i,j,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)) &
+ & -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))
+ 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.f b/src/poisson/ppm_poisson_init_predef.f
new file mode 100644
index 0000000..5baa4c8
--- /dev/null
+++ b/src/poisson/ppm_poisson_init_predef.f
@@ -0,0 +1,637 @@
+ !-------------------------------------------------------------------------
+ ! 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 __NCOM 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
+ 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
+ !@ 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) :: dx,dy,dz
+ REAL(__PREC) :: Lx2,Ly2,Lz2
+
+ 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
+
+ !-------------------------------------------------------------------------
+ ! 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)
+
+ !-------------------------------------------------------------------------
+ ! Setup mesh sizes for intermediate meshes/topologies
+ !-------------------------------------------------------------------------
+ IF (green .EQ. ppm_poisson_grn_pois_per) THEN
+ !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 !No keep first component as the size of the complex and real arrays differ - and remember to keep ndataxyc(1) etc
+ 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)
+ !Inverse of the size of the domain squared
+ Lx2 = 1.0_MK/(topology%max_physd(1)-topology%min_physd(1))**2
+ Ly2 = 1.0_MK/(topology%max_physd(2)-topology%min_physd(2))**2
+ Lz2 = 1.0_MK/(topology%max_physd(3)-topology%min_physd(3))**2
+ ELSE IF (green .EQ. ppm_poisson_grn_pois_fre) THEN !vertex
+ !Copy size of global mesh
+ ppmpoisson%nmxy (1) = mesh%nm(1)*2
+ ppmpoisson%nmxy (2) = mesh%nm(2)*2
+ ppmpoisson%nmxy (3) = mesh%nm(3)*2
+ !ppmpoisson%nmxyc(1) = (mesh%nm(1)-1)/2+1
+ ppmpoisson%nmxyc(1) = (mesh%nm(1)*2) !tmp until meshid>1 works !No keep first component as the size of the complex and real arrays differ - and remember to keep ndataxyc(1) etc
+ ppmpoisson%nmxyc(2) = mesh%nm(2)*2
+ ppmpoisson%nmxyc(3) = mesh%nm(3)*2
+ !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)*2) !tmp until meshid>1 works
+ ppmpoisson%nmz (2) = mesh%nm(2)*2
+ ppmpoisson%nmz (3) = mesh%nm(3)*2
+ !Determine the grid spacing !vertex
+ dx = (topology%max_physd(1)-topology%min_physd(1))/(mesh%nm(1)-1)
+ dy = (topology%max_physd(2)-topology%min_physd(2))/(mesh%nm(2)-1)
+ dz = (topology%max_physd(3)-topology%min_physd(3))/(mesh%nm(3)-1)
+ !maybe they should look like this: - no I think not for declaring the Greens function - after all the above dx is the value we have
+ !dx = (topology%max_physd(1)-topology%min_physd(1))/(mesh%nm(1))
+ !dy = (topology%max_physd(2)-topology%min_physd(2))/(mesh%nm(2))
+ !dz = (topology%max_physd(3)-topology%min_physd(3))/(mesh%nm(3))
+ ENDIF
+
+ !-------------------------------------------------------------------------
+ ! Create temporary derivation arrays if necessary
+ ! or spectral scaling coefficients
+ !-------------------------------------------------------------------------
+ IF (PRESENT(derive)) THEN
+ IF (( derive .EQ. ppm_poisson_drv_curl_fd2 &
+ & .OR. derive .EQ. ppm_poisson_drv_curl_fd4)) 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)))
+ !IF (( derive .EQ. ppm_poisson_drv_curl_fd2 &
+ !& .OR. derive .EQ. ppm_poisson_drv_curl_fd4) &
+ !& .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_curl_fd2 &
+ !& .OR. derive .EQ. ppm_poisson_drv_curl_fd4) &
+ !& .AND. ppmpoisson%case .EQ. ppm_poisson_grn_pois_fre) 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_curl_sp) THEN
+ ppmpoisson%normkx = 2*PI/(topology%max_physd(1)-topology%min_physd(1))
+ ppmpoisson%normky = 2*PI/(topology%max_physd(2)-topology%min_physd(2))
+ ppmpoisson%normkz = 2*PI/(topology%max_physd(3)-topology%min_physd(3))
+ ppmpoisson%derivatives = ppm_poisson_drv_curl_sp
+ 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
+
+
+
+ !-------------------------------------------------------------------------
+ ! Create new slab topology !@sofar periodic
+ !-------------------------------------------------------------------------
+ ttopoid = 0
+ tmeshid = -1
+ decomposition = ppm_param_decomp_xy_slab
+ assigning = ppm_param_assign_internal
+ IF (ppmpoisson%case .EQ. ppm_poisson_grn_pois_per) THEN
+ bcdef = ppm_param_bcdef_periodic
+ ELSE IF (ppmpoisson%case .EQ. ppm_poisson_grn_pois_fre) THEN
+ bcdef = ppm_param_bcdef_freespace
+ ENDIF
+ 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
+ IF (ppmpoisson%case .EQ. ppm_poisson_grn_pois_per) THEN
+ bcdef = ppm_param_bcdef_periodic
+ ELSE IF (ppmpoisson%case .EQ. ppm_poisson_grn_pois_fre) THEN
+ bcdef = ppm_param_bcdef_freespace
+ ENDIF
+ 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)
+
+ !-------------------------------------------------------------------------
+ ! The complex Greens function is always kept on the z-pencil topology
+ !-------------------------------------------------------------------------
+ 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
+ !!!The z-dimension is defined to full extent for the real temporary
+ !!!Greens function array
+ !!ALLOCATE(ppmpoisson%fldgrnr(&
+ !!& indl(1):indu(1),indl(2):indu(2),indl(3):((indu(3)-1)*2),&
+ !!& 1:ppmpoisson%nsublistz),stat=info)
+ ALLOCATE(ppmpoisson%fldgrnc(&
+ & 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
+ !-------------------------------------------------------------------------
+ IF (green .EQ. ppm_poisson_grn_pois_per) THEN
+ ! Scaling the spectral coefficients... one minus is due to (i*k)^2 ... how about one due to Poisson?
+ normfac = 1.0_MK/(4.0_MK*PI*PI * &
+ !and normalisation of FFTs
+ & REAL((ppmpoisson%nmz(1)-1)*(ppmpoisson%nmz(2)-1)*(ppmpoisson%nmz(3)-1),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)
+ 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...:
+ !(we subtract nm '-1' as not all points are included (periodic))
+ IF (kx .GT. (ppmpoisson%nmz(1)-1)/2) kx = kx-(ppmpoisson%nmz(1)-1)
+ IF (ky .GT. (ppmpoisson%nmz(2)-1)/2) ky = ky-(ppmpoisson%nmz(2)-1)
+ IF (kz .GT. (ppmpoisson%nmz(3)-1)/2) kz = kz-(ppmpoisson%nmz(3)-1)
+ ppmpoisson%fldgrnr(i,j,k,isub) = &
+ & normfac/(REAL(kx*kx,__PREC)*Lx2 &
+ & + REAL(ky*ky,__PREC)*Ly2 &
+ & + REAL(kz*kz,__PREC)*Lz2)
+ !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
+
+ !-------------------------------------------------------------------------
+ ! Compute real free space Greens function. Analytic
+ ! The Greens function is initialised temporarily in the real xy slabs,
+ ! then FFTed in those directions, mapped to z-pencils, FFTed in z and
+ ! finally copied to ppmpoisson%fldgrnc. The real xy slabs have already
+ ! been setup for FFTs etc so they offer a convenient container for the
+ ! FFTing the Green's function instead of setting up the whole apparatus
+ ! for this one-time affair.
+ ! These loops must run over the padded(extended) domain thus %ndataxy
+ ! \nabla \Psi = -\omega
+ ! The Greens function is 1/4piR (R is distance to origo) and includes the
+ ! minus of the RHS
+ !-------------------------------------------------------------------------
+ ELSE IF (green .EQ. ppm_poisson_grn_pois_fre) THEN
+ !-----------------------------------------------------------------------
+ ! First initialise the real Green's function
+ !@alternatively this could come from as input
+ !-----------------------------------------------------------------------
+ !@write(*,*) 'what the fuck?'
+ !there should NOT be a minus here since this Greens function takes
+ !the minus of the Poisson equation into account
+ normfac = 1.0_MK/(4.0_MK*PI* &
+ !remembering FFT normalization of ALL points: !vertex
+ & REAL((ppmpoisson%nmxy(1))*(ppmpoisson%nmxy(2))*(ppmpoisson%nmxy(3)),MK))*dx*dy*dz
+ !& REAL((ppmpoisson%nmxy(1))*(ppmpoisson%nmxy(2))*(ppmpoisson%nmxy(3)),MK)& !this is the correct normalization to bring one field back and forth.
+ !remembering FFT normalization of ALL points: !vertex
+ !!& REAL((ppmpoisson%nmxy(1))*(ppmpoisson%nmxy(2))*(ppmpoisson%nmxy(3)),MK)& !this should be correct normalization. When back and forth transforming the green's function is correct
+ !!& *REAL((ppmpoisson%nmxy(1))*(ppmpoisson%nmxy(2))*(ppmpoisson%nmxy(3))/8,MK)) !this line is probably not necessary
+ !!!& *REAL((ppmpoisson%nmxy(1))*(ppmpoisson%nmxy(2))*(ppmpoisson%nmxy(3)),MK)) !this line is probably not necessary
+ !@write(*,*) ppmpoisson%nmxy, 'johannes'
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)
+ DO j=1,ppmpoisson%ndataxy(2,isubl)
+ DO i=1,ppmpoisson%ndataxy(1,isubl)
+ kx = i-1 + (ppmpoisson%istartxy(1,isubl)-1)
+ ky = j-1 + (ppmpoisson%istartxy(2,isubl)-1)
+ kz = k-1 + (ppmpoisson%istartxy(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...:
+ !(we subtract all nmxy as all points are included in fresspace)
+ IF (kx .GT. (ppmpoisson%nmxy(1)-2)/2) kx = kx-(ppmpoisson%nmxy(1))
+ IF (ky .GT. (ppmpoisson%nmxy(2)-2)/2) ky = ky-(ppmpoisson%nmxy(2))
+ IF (kz .GT. (ppmpoisson%nmxy(3)-2)/2) kz = kz-(ppmpoisson%nmxy(3))
+ ppmpoisson%fldxyr(1,i,j,k,isub) = &
+ & normfac/(SQRT( REAL(kx*kx,__PREC)*dx*dx+ &
+ & REAL(ky*ky,__PREC)*dy*dy+ &
+ & REAL(kz*kz,__PREC)*dz*dz))
+ ppmpoisson%fldxyr(2,i,j,k,isub) = 0.0_MK
+ ppmpoisson%fldxyr(3,i,j,k,isub) = 0.0_MK
+ !Take care of singularity
+ !This is nasty as well
+ IF ((kx*kx+ky*ky+kz*kz) .EQ. 0) THEN
+ !ppmpoisson%fldxyr(1,i,j,k,isub) = 2.0_MK*normfac/(dx) - normfac/(2.0_MK*dx)
+ !ppmpoisson%fldxyr(1,i,j,k,isub) = normfac/(dx)
+ ppmpoisson%fldxyr(1,i,j,k,isub) = 0.0_MK
+ !ppmpoisson%fldxyr(1,i,j,k,isub) = 2.0_MK*normfac/(dx) - normfac/(2.0_MK*dx) + 1.0_MK
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ !-------------------------------------------------------------------------
+ ! FOURIER TRANSFORM AND MAP GALORE
+ ! @this should be updated when meshid>1 works
+ ! This part should be used both for freespace and a custom Greens function
+ !-------------------------------------------------------------------------
+ !-----------------------------------------------------------------------
+ ! Do slab FFT (XY) - using the non-reduced topology !@what does this mean
+ !-----------------------------------------------------------------------
+ CALL ppm_fft_execute_2d(ppmpoisson%topoidxy,&
+ & ppmpoisson%meshidxy, ppmpoisson%planfxy, &
+ & ppmpoisson%fldxyr, ppmpoisson%fldxyc, &
+ & info)
+ !-----------------------------------------------------------------------
+ ! 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.',isub)
+ GOTO 9999
+ ENDIF
+
+ !Push the data
+ CALL ppm_map_field_push(&
+ & ppmpoisson%topoidxy, &
+ & ppmpoisson%meshidxyc,ppmpoisson%fldxyc,__NCOM,info) !@I believe it should say meshidxy here until meshid>1 works. Then the following line will be used
+ !& 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.',isub)
+ 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.',isub)
+ 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.',isub)
+ GOTO 9999
+ ENDIF
+
+ !-----------------------------------------------------------------------
+ ! Do pencil FFT (Z)
+ !-----------------------------------------------------------------------
+ CALL ppm_fft_execute_1d(ppmpoisson%topoidz,&
+ & ppmpoisson%meshidz, ppmpoisson%planfz, &
+ & ppmpoisson%fldzc1, ppmpoisson%fldzc2, &
+ & info)
+
+ !-----------------------------------------------------------------------
+ ! Copy first component of the Fourier transformed vector to fldgrnc
+ !-----------------------------------------------------------------------
+ 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%fldgrnc(i,j,k,isub) = ppmpoisson%fldzc2(1,i,j,k,isub)!+1.0_MK
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+
+ !@Deallocate field green r NOPE!
+ 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.f b/src/poisson/ppm_poisson_solve.f
new file mode 100644
index 0000000..585fc0f
--- /dev/null
+++ b/src/poisson/ppm_poisson_solve.f
@@ -0,0 +1,837 @@
+ !-------------------------------------------------------------------------
+ ! 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) :: phix,phiy,phiz
+ 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
+!@write(*,*) 'solving .... 1', ppm_rank
+
+ !-------------------------------------------------------------------------
+ !@ Perhaps check if ghostlayer suffices for a given fd stencil
+ !-------------------------------------------------------------------------
+
+
+ !-------------------------------------------------------------------------
+ ! Perhaps allocate (and deallocate) arrays !@
+ !-------------------------------------------------------------------------
+#ifdef __VARMESH
+ !@tmp only works for one subdomain (not parallel)
+ write(*,*) 'fieldin', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)/2-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)/2-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)/2-1
+ write(*,'(A,E12.4,A,$)') ' (',fieldin(1,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+#endif
+ !-----------------------------------------------------------------------
+ ! Set the real xy slabs 0 (part of the 0 padding) for free-space
+ !@ free-space calculations and reprojection may cause problems !why?
+ !-----------------------------------------------------------------------
+ IF (presentcase .EQ. ppm_poisson_grn_pois_fre) THEN
+ ppmpoisson%fldxyr = 0.0_MK
+ ENDIF
+ !-----------------------------------------------------------------------
+ ! Map data globally to the slabs (XY)
+ ! This is where the vorticity is extended and padded with 0 for free-space
+ !-----------------------------------------------------------------------
+ !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
+!@write(*,*) 'solving .... 2', ppm_rank
+#ifdef __VARMESH
+!@tmp
+ write(*,*) 'fldxyr', 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
+#endif
+
+ !-----------------------------------------------------------------------
+ ! Do slab FFT (XY) - use the non-reduced topology !@what does this mean
+ !-----------------------------------------------------------------------
+ CALL ppm_fft_execute_2d(ppmpoisson%topoidxy,&
+ & ppmpoisson%meshidxy, ppmpoisson%planfxy, &
+ & ppmpoisson%fldxyr, ppmpoisson%fldxyc, &
+ & info)
+!@write(*,*) 'solving .... 3', ppm_rank
+
+#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
+!@write(*,*) 'solving .... 4', ppm_rank
+
+ !-----------------------------------------------------------------------
+ ! 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
+!@write(*,*) 'solving .... 5', ppm_rank
+
+ !-----------------------------------------------------------------------
+ ! Do pencil FFT (Z)
+ !-----------------------------------------------------------------------
+ CALL ppm_fft_execute_1d(ppmpoisson%topoidz,&
+ & ppmpoisson%meshidz, ppmpoisson%planfz, &
+ & ppmpoisson%fldzc1, ppmpoisson%fldzc2, &
+ & info)
+
+!@write(*,*) 'solving .... 6', ppm_rank
+#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
+!@write(*,*) 'solving .... 7', ppm_rank
+
+ !-----------------------------------------------------------------------
+ ! Apply the periodic Green's function
+ !-----------------------------------------------------------------------
+ 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
+ !-----------------------------------------------------------------------
+ ! Apply the free-space Green's function
+ !-----------------------------------------------------------------------
+ ELSE IF (presentcase .EQ. ppm_poisson_grn_pois_fre) 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%fldgrnc( i,j,k,isub)*&
+ & ppmpoisson%fldzc2(1,i,j,k,isub)
+ ppmpoisson%fldzc2(2,i,j,k,isub) = ppmpoisson%fldgrnc( i,j,k,isub)*&
+ & ppmpoisson%fldzc2(2,i,j,k,isub)
+ ppmpoisson%fldzc2(3,i,j,k,isub) = ppmpoisson%fldgrnc( i,j,k,isub)*&
+ & ppmpoisson%fldzc2(3,i,j,k,isub)
+ !!ppmpoisson%fldzc2(1,i,j,k,isub) = ppmpoisson%fldgrnc( i,j,k,isub)
+ !!ppmpoisson%fldzc2(2,i,j,k,isub) = ppmpoisson%fldgrnc( i,j,k,isub)
+ !!ppmpoisson%fldzc2(3,i,j,k,isub) = ppmpoisson%fldgrnc( i,j,k,isub)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+!@write(*,*) 'solving .... 8', ppm_rank
+ !-----------------------------------------------------------------------
+ ! Spectral derivatives
+ ! normkx, etc contains 2pi/Lx
+ !-----------------------------------------------------------------------
+ IF (ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_sp .AND.&
+ & (presentcase .EQ. ppm_poisson_grn_pois_per .OR. &
+ presentcase .EQ. ppm_poisson_grn_pois_fre)) THEN
+ write(*,*) 'curl fft style! not working yet though!'
+ !remembering to normalize the FFT
+ normfac = 1.0_MK/ REAL((ppmpoisson%nmz(1)-1)* & !vertex
+ & (ppmpoisson%nmz(2)-1)* &
+ & (ppmpoisson%nmz(3)-1),MK)
+ 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 = CMPLX(0.0_MK,gk*ppmpoisson%normkz,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 = CMPLX(0.0_MK,gj*ppmpoisson%normky,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 = CMPLX(0.0_MK,gi*ppmpoisson%normkx,MK)
+
+ phix = ppmpoisson%fldzc2(1,i,j,k,isub)
+ phiy = ppmpoisson%fldzc2(2,i,j,k,isub)
+ phiz = ppmpoisson%fldzc2(3,i,j,k,isub)
+
+ !maybe normfac on kx,ky,kz?
+ ppmpoisson%fldzc2(1,i,j,k,isub) = normfac*(ky*phiz-kz*phiy)
+ ppmpoisson%fldzc2(2,i,j,k,isub) = normfac*(kz*phix-kx*phiz)
+ ppmpoisson%fldzc2(3,i,j,k,isub) = normfac*(kx*phiy-ky*phix)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+
+ !-----------------------------------------------------------------------
+ ! Vorticity re-projection
+ !@ has the domain length really been included in these wave numbers
+ !-----------------------------------------------------------------------
+ ELSE IF (presentcase .EQ. ppm_poisson_grn_reprojec) THEN
+ !remembering to normalize the FFT
+ normfac = 1.0_MK/ REAL((ppmpoisson%nmz(1)-1)* & !vertex
+ & (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
+
+ 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
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+
+!@write(*,*) 'solving .... 9', ppm_rank
+ !-----------------------------------------------------------------------
+ ! IFFT pencil (Z)
+ !-----------------------------------------------------------------------
+ CALL ppm_fft_execute_1d(ppmpoisson%topoidz,&
+ & ppmpoisson%meshidz, ppmpoisson%planbz, &
+ & ppmpoisson%fldzc2, ppmpoisson%fldzc1, &
+ & info)
+
+!@write(*,*) 'solving .... 11', ppm_rank
+#ifdef __NOPE
+ if (ppm_rank .EQ. trank) THEN
+ write(*,*)
+ write(*,*)
+ 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
+
+ !-----------------------------------------------------------------------
+ ! 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
+
+ !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(*,*) 'solving .... 12', ppm_rank
+ !-----------------------------------------------------------------------
+ ! IFFT (XY) use the non-reduced topology
+ !-----------------------------------------------------------------------
+ CALL ppm_fft_execute_2d(ppmpoisson%topoidxy,&
+ & ppmpoisson%meshidxy, ppmpoisson%planbxy, &
+ & ppmpoisson%fldxyc, ppmpoisson%fldxyr, &
+ & info)
+
+!@write(*,*) 'solving .... 13', ppm_rank
+#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
+
+#ifdef __VARMESH
+!@tmp
+ write(*,*) 'fldxyr2', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)
+ DO j=1,ppmpoisson%ndataxy(2,isubl)
+ write(*,'(A,E12.4,A,$)') ' (',ppmpoisson%fldxyr(1,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+#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
+!@write(*,*) 'solving .... 13a', ppm_rank
+
+ !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
+!@write(*,*) 'solving .... 13b', ppm_rank, ppmpoisson%topoidxy, topoid, ppmpoisson%meshidxy, meshid
+
+ !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
+
+!@write(*,*) 'solving .... 14', ppm_rank
+ !-------------------------------------------------------------------------
+ ! 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 .OR. &
+ ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_fd4) .AND. &
+ & (presentcase .EQ. ppm_poisson_grn_pois_per .OR. & !@these may be unnecessary - perhaps just the derive value. Or maybe not: in case of vorticity reprojection we could get lost
+ presentcase .EQ. ppm_poisson_grn_pois_fre)) THEN
+!@write(*,*) 'solving .... 14a', ppm_rank
+ CALL ppm_map_field_pop(&
+ & topoid, &
+ & meshid,ppmpoisson%drv_vr, &
+ & __NCOM,gstw,info) !!! gst
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
+ GOTO 9999
+ ENDIF
+!@write(*,*) 'solving .... 14b', ppm_rank
+ !-------------------------------------------------------------------------
+ ! Ghost the temporary array for derivatives (drv_vr)
+ !-------------------------------------------------------------------------
+!@write(*,*) 'solving .... 14c', ppm_rank
+ CALL ppm_map_field_ghost_get(topoid,meshid,gstw,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to initialise ghosts.',info2)
+ GOTO 9999
+ ENDIF
+!@write(*,*) 'solving .... 14d', ppm_rank
+ CALL ppm_map_field_push(topoid,meshid,ppmpoisson%drv_vr,__NCOM,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to push ghosts.',info2)
+ GOTO 9999
+ ENDIF
+!@write(*,*) 'solving .... 14e', ppm_rank
+ CALL ppm_map_field_send(info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to send ghosts.',info2)
+ GOTO 9999
+ ENDIF
+!@write(*,*) 'solving .... 14f', ppm_rank
+ CALL ppm_map_field_pop(topoid,meshid,ppmpoisson%drv_vr,__NCOM,gstw,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop ghosts.',info2)
+ GOTO 9999
+ ENDIF
+
+ ELSE
+!@write(*,*) 'solving .... 14g', ppm_rank
+ CALL ppm_map_field_pop(&
+ & topoid, &
+ & meshid,fieldout, &
+ & __NCOM,gstw,info) !!! gst
+!@write(*,*) 'solving .... 14h', ppm_rank
+ ENDIF
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
+ GOTO 9999
+ ENDIF
+
+!@write(*,*) 'solving .... 15', ppm_rank
+#ifdef __VARMESH
+!@tmp only works for one subdomain (not parallel)
+ write(*,*) 'fieldout', ppm_rank
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ DO k=1,ppmpoisson%ndataxy(3,isubl)/2-1
+ write(*,*) 'z',k
+ DO i=1,ppmpoisson%ndataxy(1,isubl)/2-1
+ DO j=1,ppmpoisson%ndataxy(2,isubl)/2-1
+ write(*,'(A,E12.4,A,$)') ' (',fieldout(1,i,j,k,isub),')'
+ END DO
+ write(*,*)
+ END DO
+ END DO
+ END DO
+#endif
+
+ !-------------------------------------------------------------------------
+ !@To come: treat ghost layer to make FD stencils work
+ !-------------------------------------------------------------------------
+ IF (ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_fd2 .AND.&
+ & (presentcase .EQ. ppm_poisson_grn_pois_fre)) THEN
+ CALL ppm_poisson_extrapolateghost(topoid,meshid,ppmpoisson%drv_vr,&
+ & 2,4,gstw,info)
+ ENDIF
+ IF (ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_fd4 .AND.&
+ & (presentcase .EQ. ppm_poisson_grn_pois_fre)) THEN
+ CALL ppm_poisson_extrapolateghost(topoid,meshid,ppmpoisson%drv_vr,&
+ & 2,4,gstw,info)
+ ENDIF
+
+!@write(*,*) 'solving .... 16', ppm_rank
+ !-------------------------------------------------------------------------
+ ! 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 .OR. & !@these may be unnecessary - perhaps just the derive value
+ presentcase .EQ. ppm_poisson_grn_pois_fre)) THEN
+ CALL ppm_poisson_fd(topoid,meshid,ppmpoisson%drv_vr,fieldout,&
+ & ppm_poisson_drv_curl_fd2,info)
+ ENDIF
+ IF (ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_fd4 .AND.&
+ & (presentcase .EQ. ppm_poisson_grn_pois_per .OR. & !@these may be unnecessary - perhaps just the derive value
+ presentcase .EQ. ppm_poisson_grn_pois_fre)) THEN
+ CALL ppm_poisson_fd(topoid,meshid,ppmpoisson%drv_vr,fieldout,&
+ & ppm_poisson_drv_curl_fd4,info)
+ ENDIF
+
+!@write(*,*) 'solving .... 17', ppm_rank
+ !-------------------------------------------------------------------------
+ ! Finally ghost the velocity/stream function field before returning it
+ ! Also extrapolate if freespace
+ !-------------------------------------------------------------------------
+ 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)
+ IF (presentcase .EQ. ppm_poisson_grn_pois_fre) THEN
+ CALL ppm_poisson_extrapolateghost(topoid,meshid,fieldout,&
+ & 2,4,gstw,info)
+ ENDIF
+
+!@write(*,*) 'solving .... 18', ppm_rank
+ !-------------------------------------------------------------------------
+ ! Perhaps allocate (and deallocate) arrays !@
+ !-------------------------------------------------------------------------
+
+ !-------------------------------------------------------------------------
+ ! Return
+ !-------------------------------------------------------------------------
+ 9999 CONTINUE
+ CALL substop('ppm_poisson_solve',t0,info)
+ RETURN
+
+ END SUBROUTINE __ROUTINE
+
--
GitLab