Skip to content
Snippets Groups Projects
Commit 65c32a78 authored by Omar Awile's avatar Omar Awile
Browse files

Cleanup 

Cleaned up all modules and subroutines that will not be immediately
ported to the new PPM core code base. Renamed and combined some ODE
modules.
parent c54e8f1f
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 8366 deletions
src/fft/ppm_fft_exec_3d_vec_bc2r_xy.f deleted 100644 → 0
+ 0
− 121
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_fft_exec_3d_vec_bc2r_xy
!-------------------------------------------------------------------------
! Copyright (c) 2012 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
#undef __ROUTINE
#undef __PREC
src/fft/ppm_fft_exec_3d_vec_c2c_z.f deleted 100644 → 0
+ 0
− 126
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_fft_exec_3d_vec_c2c_z
!-------------------------------------------------------------------------
! Copyright (c) 2012 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
#undef __ROUTINE
#undef __PREC
src/fft/ppm_fft_exec_3d_vec_fr2c_xy.f deleted 100644 → 0
+ 0
− 105
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_fft_exec_3d_vec_fr2c_xy
!-------------------------------------------------------------------------
! Copyright (c) 2012 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
#undef __ROUTINE
#undef __PREC
src/fft/ppm_fft_normalize_c.f deleted 100644 → 0
+ 0
− 129
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_fft_normalize_c
!-------------------------------------------------------------------------
! Copyright (c) 2012 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
#undef __ROUTINE
#undef __PREC
src/fft/ppm_fft_normalize_r.f deleted 100644 → 0
+ 0
− 130
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_fft_normalize_r
!-------------------------------------------------------------------------
! Copyright (c) 2012 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
#undef __ROUTINE
#undef __PREC
src/fft/ppm_fft_plan_3d_vec_bc2c_z.f deleted 100644 → 0
+ 0
− 151
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_fft_plan_3d_vec_bc2c_z
!-------------------------------------------------------------------------
! Copyright (c) 2012 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%nm(3)-1
ELSE
ppmplan%nx(1,isub) = mesh%nm(3)
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
#undef __ROUTINE
#undef __PREC
src/fft/ppm_fft_plan_3d_vec_bc2r_xy.f deleted 100644 → 0
+ 0
− 156
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_fft_plan_3d_vec_bc2r_xy
!-------------------------------------------------------------------------
! Copyright (c) 2012 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
!the number of points in each direction of the piece to be transformed
!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%nm(1)-1
ppmplan%nx(2,isub) = mesh%nm(2)-1
ELSE
ppmplan%nx(1,isub) = mesh%nm(1)
ppmplan%nx(2,isub) = mesh%nm(2)
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
#undef __ROUTINE
#undef __PREC
src/fft/ppm_fft_plan_3d_vec_fc2c_z.f deleted 100644 → 0
+ 0
− 151
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_fft_plan_3d_vec_fc2c_z
!-------------------------------------------------------------------------
! Copyright (c) 2012 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%nm(3)-1
ELSE
ppmplan%nx(1,isub) = mesh%nm(3)
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
#undef __ROUTINE
#undef __PREC
src/fft/ppm_fft_plan_3d_vec_fr2c_xy.f deleted 100644 → 0
+ 0
− 155
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_fft_plan_3d_vec_fr2c_xy
!-------------------------------------------------------------------------
! Copyright (c) 2012 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%nm(1)-1
ppmplan%nx(2,isub) = mesh%nm(2)-1
ELSE
ppmplan%nx(1,isub) = mesh%nm(1)
ppmplan%nx(2,isub) = mesh%nm(2)
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
#undef __ROUTINE
#undef __PREC
src/mg/ppm_mg_alloc_bc.f deleted 100644 → 0
+ 0
− 396
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_mg_alloc_bc
!-------------------------------------------------------------------------
!
! Purpose : Does the (re)allocation of arrays of type
! bc_value. It offers the same allocation
! types as ppm_alloc for regular arrays.
!
! Input : lda(:) (I) number of subdomains +
! numbers of levels
! iopt (I) alloc action. One of:
! ppm_param_alloc_fit_preserve
! ppm_param_alloc_fit
! ppm_param_alloc_grow_preserve
! ppm_param_alloc_grow
! ppm_param_dealloc
!
! Input/output : field (T) array of TYPE(bc_value)
! which is to be (re)allocated.
!
! Output : info (I) Return status. 0 if everything OK.
!
! Remarks :
!
! References :
!
! Revisions :
!-------------------------------------------------------------------------
! $Log: ppm_mg_alloc_bc.f,v $
! Revision 1.1.1.1 2007/07/13 10:18:56 ivos
! CBL version of the PPM library
!
! Revision 1.7 2004/10/01 16:33:38 ivos
! cosmetics.
!
! Revision 1.6 2004/10/01 16:09:10 ivos
! Replaced REAL(ppm_kind_double) :: t0 with REAL(MK) t0.
!
! Revision 1.5 2004/09/22 18:25:08 kotsalie
! MG new version
!
! Revision 1.4 2004/07/26 13:49:18 ivos
! Removed Routines sections from the header comment.
!
! Revision 1.3 2004/07/26 11:59:39 ivos
! Fixes to make it compile.
!
! Revision 1.2 2004/07/26 07:46:37 ivos
! Changed to use single-interface modules. Updated all USE statements.
!
! Revision 1.1 2004/06/29 14:43:14 kotsalie
! Needed for my type allocation
!
!-------------------------------------------------------------------------
! Copyright (c) 2012 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
! Center for Fluid Dynamics (DTU)
!
!
! This file is part of the Parallel Particle Mesh Library (PPM).
!
! PPM is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License
! as published by the Free Software Foundation, either
! version 3 of the License, or (at your option) any later
! version.
!
! PPM is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! and the GNU Lesser General Public License along with PPM. If not,
! see <http://www.gnu.org/licenses/>.
!
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_alloc_bc_2d_sca_s(field,lda,iopt,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_alloc_bc_2d_sca_d(field,lda,iopt,info)
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_alloc_bc_3d_sca_s(field,lda,iopt,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_alloc_bc_3d_sca_d(field,lda,iopt,info)
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_alloc_bc_2d_vec_s(field,lda,iopt,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_alloc_bc_2d_vec_d(field,lda,iopt,info)
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_alloc_bc_3d_vec_s(field,lda,iopt,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_alloc_bc_3d_vec_d(field,lda,iopt,info)
#endif
#endif
#endif
!-------------------------------------------------------------------------
! Includes
!-------------------------------------------------------------------------
#include "ppm_define.h"
!-------------------------------------------------------------------------
! Module
!-------------------------------------------------------------------------
USE ppm_module_data
USE ppm_module_data_mg
USE ppm_module_substart
USE ppm_module_substop
USE ppm_module_error
USE ppm_module_alloc
USE ppm_module_write
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_single
#else
INTEGER, PARAMETER :: MK = ppm_kind_double
#endif
!-------------------------------------------------------------------------
! Arguments
!-------------------------------------------------------------------------
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND ==__SINGLE_PRECISION
TYPE(bc_value_2d_sca_s), DIMENSION(:), POINTER :: field
#elif __KIND == __DOUBLE_PRECISION
TYPE(bc_value_2d_sca_d), DIMENSION(:), POINTER :: field
#endif
#elif __MESH_DIM == __3D
#if __KIND ==__SINGLE_PRECISION
TYPE(bc_value_3d_sca_s), DIMENSION(:), POINTER :: field
#elif __KIND == __DOUBLE_PRECISION
TYPE(bc_value_3d_sca_d), DIMENSION(:), POINTER :: field
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND ==__SINGLE_PRECISION
TYPE(bc_value_2d_vec_s), DIMENSION(:), POINTER :: field
#elif __KIND == __DOUBLE_PRECISION
TYPE(bc_value_2d_vec_d), DIMENSION(:), POINTER :: field
#endif
#elif __MESH_DIM == __3D
#if __KIND ==__SINGLE_PRECISION
TYPE(bc_value_3d_vec_s), DIMENSION(:), POINTER :: field
#elif __KIND == __DOUBLE_PRECISION
TYPE(bc_value_3d_vec_d), DIMENSION(:), POINTER :: field
#endif
#endif
#endif
INTEGER, DIMENSION(: ), INTENT(IN ) :: lda
INTEGER, INTENT(IN ) :: iopt
INTEGER, INTENT( OUT) :: info
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
INTEGER :: i
INTEGER, DIMENSION(1) :: ldc
REAL(MK) :: t0
LOGICAL :: lcopy,lalloc,lrealloc,ldealloc
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND ==__SINGLE_PRECISION
TYPE(bc_value_2d_sca_s), DIMENSION(:), POINTER :: work_field
#elif __KIND == __DOUBLE_PRECISION
TYPE(bc_value_2d_sca_d), DIMENSION(:), POINTER :: work_field
#endif
#elif __MESH_DIM == __3D
#if __KIND ==__SINGLE_PRECISION
TYPE(bc_value_3d_sca_s), DIMENSION(:), POINTER :: work_field
#elif __KIND == __DOUBLE_PRECISION
TYPE(bc_value_3d_sca_d), DIMENSION(:), POINTER :: work_field
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND ==__SINGLE_PRECISION
TYPE(bc_value_2d_vec_s), DIMENSION(:), POINTER :: work_field
#elif __KIND == __DOUBLE_PRECISION
TYPE(bc_value_2d_vec_d), DIMENSION(:), POINTER :: work_field
#endif
#elif __MESH_DIM == __3D
#if __KIND ==__SINGLE_PRECISION
TYPE(bc_value_3d_vec_s), DIMENSION(:), POINTER :: work_field
#elif __KIND == __DOUBLE_PRECISION
TYPE(bc_value_3d_vec_d), DIMENSION(:), POINTER :: work_field
#endif
#endif
#endif
!-------------------------------------------------------------------------
! Externals
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Initialise
!-------------------------------------------------------------------------
CALL substart('ppm_mg_alloc_bc',t0,info)
!-------------------------------------------------------------------------
! Check arguments
!-------------------------------------------------------------------------
IF (ppm_debug .GT. 0 .AND. iopt .NE. ppm_param_dealloc) THEN
IF (SIZE(lda,1) .LT. 1) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_mg_alloc_field', &
& 'lda must be at least of length 1',__LINE__,info)
GOTO 9999
ENDIF
IF (lda(1) .LT. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_mg_alloc_field', &
& 'lda(1) must be >= 0',__LINE__,info)
GOTO 9999
ENDIF
ENDIF
!-------------------------------------------------------------------------
! Check allocation type
!-------------------------------------------------------------------------
lcopy = .FALSE.
lalloc = .FALSE.
lrealloc = .FALSE.
ldealloc = .FALSE.
IF (iopt .EQ. ppm_param_alloc_fit_preserve) THEN
!---------------------------------------------------------------------
! Fit memory and preserve the present contents
!---------------------------------------------------------------------
IF (ASSOCIATED(field)) THEN
ldc(1) = SIZE(field,1)
IF (ldc(1) .NE. lda(1)) THEN
lalloc = .TRUE.
lrealloc = .TRUE.
lcopy = .TRUE.
ENDIF
ELSE
lalloc = .TRUE.
ENDIF
ELSEIF (iopt .EQ. ppm_param_alloc_fit) THEN
!---------------------------------------------------------------------
! Fit memory and discard the present contents
!---------------------------------------------------------------------
IF (ASSOCIATED(field)) THEN
ldc(1) = SIZE(field,1)
IF (ldc(1) .NE. lda(1)) THEN
lalloc = .TRUE.
lrealloc = .TRUE.
ldealloc = .TRUE.
ENDIF
ELSE
lalloc = .TRUE.
ENDIF
ELSEIF (iopt .EQ. ppm_param_alloc_grow_preserve) THEN
!---------------------------------------------------------------------
! Fit memory and preserve the present contents
!---------------------------------------------------------------------
IF (ASSOCIATED(field)) THEN
ldc(1) = SIZE(field,1)
IF (ldc(1) .LT. lda(1)) THEN
lalloc = .TRUE.
lrealloc = .TRUE.
lcopy = .TRUE.
ENDIF
ELSE
lalloc = .TRUE.
ENDIF
ELSEIF (iopt .EQ. ppm_param_alloc_grow) THEN
!---------------------------------------------------------------------
! Fit memory and discard the present contents
!---------------------------------------------------------------------
IF (ASSOCIATED(field)) THEN
ldc(1) = SIZE(field,1)
IF (ldc(1) .LT. lda(1)) THEN
lalloc = .TRUE.
lrealloc = .TRUE.
ldealloc = .TRUE.
ENDIF
ELSE
lalloc = .TRUE.
ENDIF
ELSEIF (iopt .EQ. ppm_param_dealloc) THEN
!---------------------------------------------------------------------
! Deallocate
!---------------------------------------------------------------------
IF (ASSOCIATED(field)) THEN
ldc(1) = SIZE(field,1)
lrealloc = .TRUE.
ldealloc = .TRUE.
ENDIF
ENDIF
!-------------------------------------------------------------------------
! Perform the actual alloc action
!-------------------------------------------------------------------------
IF (lalloc) THEN
!---------------------------------------------------------------------
! Allocate new array with new size and nullify all members
!---------------------------------------------------------------------
ALLOCATE(work_field(lda(1)),STAT=info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_mg_alloc_field', &
& 'new mgfield WORK_MESH',__LINE__,info)
GOTO 9999
ENDIF
DO i=1,lda(1)
NULLIFY(work_field(i)%pbcvalue)
ENDDO
ENDIF
IF (lcopy) THEN
!---------------------------------------------------------------------
! Save the old contents
!---------------------------------------------------------------------
DO i=1,MIN(ldc(1),lda(1))
work_field(i)%pbcvalue => field(i)%pbcvalue
ENDDO
ENDIF
IF (ldealloc) THEN
!---------------------------------------------------------------------
! Deallocate the old contents
!---------------------------------------------------------------------
DO i=1,ldc(1)
IF (ASSOCIATED(field(i)%pbcvalue)) THEN
DEALLOCATE(field(i)%pbcvalue,STAT=info)
IF (info .NE. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_dealloc,'ppm_mg_alloc_bc',&
& 'field correction FIELD%PBCVALUE',__LINE__,info)
ENDIF
NULLIFY(field(i)%pbcvalue)
ENDIF
ENDDO
ENDIF
IF (lrealloc) THEN
!---------------------------------------------------------------------
! Deallocate old pointer array
!---------------------------------------------------------------------
DEALLOCATE(field,STAT=info)
IF (info .NE. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_dealloc,'ppm_mg_alloc_bc', &
& 'mgfield data FIELD',__LINE__,info)
ENDIF
NULLIFY(field)
ENDIF
IF (lalloc) THEN
!---------------------------------------------------------------------
! Point result to new array
!---------------------------------------------------------------------
field => work_field
ENDIF
!-------------------------------------------------------------------------
! Return
!-------------------------------------------------------------------------
9999 CONTINUE
CALL substop('ppm_mg_alloc_bc',t0,info)
RETURN
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_alloc_bc_2d_sca_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_alloc_bc_2d_sca_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_alloc_bc_3d_sca_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_alloc_bc_3d_sca_d
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_alloc_bc_2d_vec_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_alloc_bc_2d_vec_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_alloc_bc_3d_vec_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_alloc_bc_3d_vec_d
#endif
#endif
#endif
src/mg/ppm_mg_alloc_field.f deleted 100644 → 0
+ 0
− 428
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_mg_alloc_field
!-------------------------------------------------------------------------
! This file is part of the Parallel Particle Mesh Library (PPM).
!
! PPM is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License
! as published by the Free Software Foundation, either
! version 3 of the License, or (at your option) any later
! version.
!
! PPM is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! and the GNU Lesser General Public License along with PPM. If not,
! see <http://www.gnu.org/licenses/>.
!
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
!
! Purpose : Does the (re)allocation of arrays of type
! mg_field. It offers the same allocation
! types as ppm_alloc for regular arrays.
!
! Input : lda(:) (I) number of subdomains +
! numbers of levels
! iopt (I) alloc action. One of:
! ppm_param_alloc_fit_preserve
! ppm_param_alloc_fit
! ppm_param_alloc_grow_preserve
! ppm_param_alloc_grow
! ppm_param_dealloc
!
! Input/output : field (T) array of TYPE(mg_field)
! which is to be (re)allocated.
!
! Output : info (I) Return status. 0 if everything OK.
!
! Remarks :
!
! References :
!
! Revisions :
!
!
!-------------------------------------------------------------------------
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_alloc_field_2d_sca_s(field,lda,iopt,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_alloc_field_2d_sca_d(field,lda,iopt,info)
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_alloc_field_3d_sca_s(field,lda,iopt,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_alloc_field_3d_sca_d(field,lda,iopt,info)
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_alloc_field_2d_vec_s(field,lda,iopt,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_alloc_field_2d_vec_d(field,lda,iopt,info)
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_alloc_field_3d_vec_s(field,lda,iopt,info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_alloc_field_3d_vec_d(field,lda,iopt,info)
#endif
#endif
#endif
!-------------------------------------------------------------------------
! Includes
!-------------------------------------------------------------------------
#include "ppm_define.h"
!-------------------------------------------------------------------------
! Module
!-------------------------------------------------------------------------
USE ppm_module_data
USE ppm_module_data_mg
USE ppm_module_substart
USE ppm_module_substop
USE ppm_module_error
USE ppm_module_write
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_single
#else
INTEGER, PARAMETER :: MK = ppm_kind_double
#endif
!-------------------------------------------------------------------------
! Arguments
!-------------------------------------------------------------------------
INTEGER , DIMENSION(: ), INTENT(IN ) :: lda
INTEGER , INTENT(IN ) :: iopt
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND ==__SINGLE_PRECISION
TYPE(mg_field_2d_sca_s), DIMENSION(:,:), POINTER :: field
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_2d_sca_d), DIMENSION(:,:), POINTER :: field
#endif
#elif __MESH_DIM == __3D
#if __KIND ==__SINGLE_PRECISION
TYPE(mg_field_3d_sca_s), DIMENSION(:,:), POINTER :: field
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_3d_sca_d),DIMENSION(:,:), POINTER :: field
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND ==__SINGLE_PRECISION
TYPE(mg_field_2d_vec_s), DIMENSION(:,:), POINTER :: field
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_2d_vec_d), DIMENSION(:,:), POINTER :: field
#endif
#elif __MESH_DIM == __3D
#if __KIND ==__SINGLE_PRECISION
TYPE(mg_field_3d_vec_s), DIMENSION(:,:), POINTER :: field
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_3d_vec_d),DIMENSION(:,:), POINTER :: field
#endif
#endif
#endif
INTEGER , INTENT( OUT) :: info
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
INTEGER :: i,j
INTEGER, DIMENSION(2) :: ldc
REAL(MK) :: t0
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND ==__SINGLE_PRECISION
TYPE(mg_field_2d_sca_s), DIMENSION(:,:), POINTER :: work_field
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_2d_sca_d), DIMENSION(:,:), POINTER :: work_field
#endif
#elif __MESH_DIM == __3D
#if __KIND ==__SINGLE_PRECISION
TYPE(mg_field_3d_sca_s), DIMENSION(:,:), POINTER :: work_field
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_3d_sca_d), DIMENSION(:,:), POINTER :: work_field
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND ==__SINGLE_PRECISION
TYPE(mg_field_2d_vec_s), DIMENSION(:,:), POINTER :: work_field
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_2d_vec_d), DIMENSION(:,:), POINTER :: work_field
#endif
#elif __MESH_DIM == __3D
#if __KIND ==__SINGLE_PRECISION
TYPE(mg_field_3d_vec_s), DIMENSION(:,:), POINTER :: work_field
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_3d_vec_d), DIMENSION(:,:), POINTER :: work_field
#endif
#endif
#endif
LOGICAL :: lcopy,lalloc,lrealloc,ldealloc
!-------------------------------------------------------------------------
! Externals
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Initialise
!-------------------------------------------------------------------------
CALL substart('ppm_mg_alloc_field',t0,info)
!-------------------------------------------------------------------------
! Check arguments
!-------------------------------------------------------------------------
IF (ppm_debug .GT. 0 .AND. iopt .NE. ppm_param_dealloc) THEN
IF (SIZE(lda,1) .LT. 2) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_mg_alloc_field', &
& 'lda must be at least of length 2',__LINE__,info)
GOTO 9999
ENDIF
IF (lda(1) .LT. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_mg_alloc_field', &
& 'lda(1) must be >= 0',__LINE__,info)
GOTO 9999
ENDIF
IF (lda(2) .LT. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_argument,'ppm_mg_alloc_field', &
& 'lda(2) must be >= 0',__LINE__,info)
GOTO 9999
ENDIF
ENDIF
!-------------------------------------------------------------------------
! Check allocation type
!-------------------------------------------------------------------------
lcopy = .FALSE.
lalloc = .FALSE.
lrealloc = .FALSE.
ldealloc = .FALSE.
IF (iopt .EQ. ppm_param_alloc_fit_preserve) THEN
!---------------------------------------------------------------------
! Fit memory and preserve the present contents
!---------------------------------------------------------------------
IF (ASSOCIATED(field)) THEN
ldc(1) = SIZE(field,1)
ldc(2) = SIZE(field,2)
IF ((ldc(1) .NE. lda(1)) .OR. (ldc(2) .NE. lda(2))) THEN
lalloc = .TRUE.
lrealloc = .TRUE.
lcopy = .TRUE.
ENDIF
ELSE
lalloc = .TRUE.
ENDIF
ELSEIF (iopt .EQ. ppm_param_alloc_fit) THEN
!---------------------------------------------------------------------
! Fit memory and discard the present contents
!---------------------------------------------------------------------
IF (ASSOCIATED(field)) THEN
ldc(1) = SIZE(field,1)
ldc(2) = SIZE(field,2)
IF ((ldc(1) .NE. lda(1)) .OR. (ldc(2) .NE. lda(2))) THEN
lalloc = .TRUE.
lrealloc = .TRUE.
ldealloc = .TRUE.
ENDIF
ELSE
lalloc = .TRUE.
ENDIF
ELSEIF (iopt .EQ. ppm_param_alloc_grow_preserve) THEN
!---------------------------------------------------------------------
! Fit memory and preserve the present contents
!---------------------------------------------------------------------
IF (ASSOCIATED(field)) THEN
ldc(1) = SIZE(field,1)
ldc(2) = SIZE(field,2)
IF ((ldc(1) .LT. lda(1)) .OR. (ldc(2) .LT. lda(2))) THEN
lalloc = .TRUE.
lrealloc = .TRUE.
lcopy = .TRUE.
ENDIF
ELSE
lalloc = .TRUE.
ENDIF
ELSEIF (iopt .EQ. ppm_param_alloc_grow) THEN
!---------------------------------------------------------------------
! Fit memory and discard the present contents
!---------------------------------------------------------------------
IF (ASSOCIATED(field)) THEN
ldc(1) = SIZE(field,1)
ldc(2) = SIZE(field,2)
IF ((ldc(1) .LT. lda(1)) .OR. (ldc(2) .LT. lda(2))) THEN
lalloc = .TRUE.
lrealloc = .TRUE.
ldealloc = .TRUE.
ENDIF
ELSE
lalloc = .TRUE.
ENDIF
ELSEIF (iopt .EQ. ppm_param_dealloc) THEN
!---------------------------------------------------------------------
! Deallocate
!---------------------------------------------------------------------
IF (ASSOCIATED(field)) THEN
ldc(1) = SIZE(field,1)
ldc(2) = SIZE(field,2)
lrealloc = .TRUE.
ldealloc = .TRUE.
ENDIF
ENDIF
!-------------------------------------------------------------------------
! Perform the actual alloc action
!-------------------------------------------------------------------------
IF (lalloc) THEN
!---------------------------------------------------------------------
! Allocate new array with new size and nullify all members
!---------------------------------------------------------------------
ALLOCATE(work_field(lda(1),lda(2)),STAT=info)
IF (info .NE. 0) THEN
info = ppm_error_fatal
CALL ppm_error(ppm_err_alloc,'ppm_mg_alloc_field', &
& 'new mgfield WORK_MESH',__LINE__,info)
GOTO 9999
ENDIF
DO j=1,lda(2)
DO i=1,lda(1)
NULLIFY(work_field(i,j)%uc)
NULLIFY(work_field(i,j)%fc)
NULLIFY(work_field(i,j)%err)
NULLIFY(work_field(i,j)%bcvalue)
ENDDO
ENDDO
ENDIF
IF (lcopy) THEN
!---------------------------------------------------------------------
! Save the old contents
!---------------------------------------------------------------------
DO j=1,MIN(ldc(2),lda(2))
DO i=1,MIN(ldc(1),lda(1))
work_field(i,j)%uc => field(i,j)%uc
work_field(i,j)%fc => field(i,j)%fc
work_field(i,j)%err => field(i,j)%err
work_field(i,j)%bcvalue => field(i,j)%bcvalue
ENDDO
ENDDO
ENDIF
IF (ldealloc) THEN
!---------------------------------------------------------------------
! Deallocate the old contents
!---------------------------------------------------------------------
DO j=1,ldc(2)
DO i=1,ldc(1)
IF (ASSOCIATED(field(i,j)%uc)) THEN
DEALLOCATE(field(i,j)%uc,STAT=info)
IF (info .NE. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_dealloc,'ppm_mg_alloc_field',&
& 'field correction FIELD%UC',__LINE__,info)
ENDIF
NULLIFY(field(i,j)%uc)
ENDIF
IF (ASSOCIATED(field(i,j)%fc)) THEN
DEALLOCATE(field(i,j)%fc,STAT=info)
IF (info .NE. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_dealloc,'ppm_mg_alloc_field',&
& 'error restriction FIELD%FC',__LINE__,info)
ENDIF
NULLIFY(field(i,j)%fc)
ENDIF
IF (ASSOCIATED(field(i,j)%err)) THEN
DEALLOCATE(field(i,j)%err,STAT=info)
IF (info .NE. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_dealloc,'ppm_mg_alloc_field',&
& 'residual FIELD%ERR',__LINE__,info)
ENDIF
NULLIFY(field(i,j)%err)
ENDIF
IF (ASSOCIATED(field(i,j)%bcvalue)) THEN
DEALLOCATE(field(i,j)%bcvalue,STAT=info)
IF (info .NE. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_dealloc,'ppm_mg_alloc_field',&
& 'boundary values FIELD%BCVALUE',__LINE__,info)
ENDIF
NULLIFY(field(i,j)%bcvalue)
ENDIF
ENDDO
ENDDO
ENDIF
IF (lrealloc) THEN
!---------------------------------------------------------------------
! Deallocate old pointer array
!---------------------------------------------------------------------
DEALLOCATE(field,STAT=info)
IF (info .NE. 0) THEN
info = ppm_error_error
CALL ppm_error(ppm_err_dealloc,'ppm_mg_alloc_field', &
& 'mgfield data FIELD',__LINE__,info)
ENDIF
NULLIFY(field)
ENDIF
IF (lalloc) THEN
!---------------------------------------------------------------------
! Point result to new array
!---------------------------------------------------------------------
field => work_field
ENDIF
!-------------------------------------------------------------------------
! Return
!-------------------------------------------------------------------------
9999 CONTINUE
CALL substop('ppm_mg_alloc_field',t0,info)
RETURN
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_alloc_field_2d_sca_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_alloc_field_2d_sca_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_alloc_field_3d_sca_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_alloc_field_3d_sca_d
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_alloc_field_2d_vec_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_alloc_field_2d_vec_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_alloc_field_3d_vec_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_alloc_field_3d_vec_d
#endif
#endif
#endif
src/mg/ppm_mg_core.f deleted 100644 → 0
+ 0
− 592
View file @ c54e8f1f
This diff is collapsed.
src/mg/ppm_mg_finalize.f deleted 100644 → 0
+ 0
− 284
View file @ c54e8f1f
!-------------------------------------------------------------------------
! Subroutine : ppm_mg_finalize
!-------------------------------------------------------------------------
!
! Purpose : This routine deallocates all the arrays
!
!
! Input :
!
! Input/output :
!
! Output : info (I) 0 on success.
!
! Remarks :
!
! References :
!
! Revisions :
!-------------------------------------------------------------------------
! $Log: ppm_mg_finalize.f,v $
! Revision 1.1.1.1 2007/07/13 10:18:56 ivos
! CBL version of the PPM library
!
! Revision 1.5 2006/07/21 11:30:55 kotsalie
! FRIDAY
!
! Revision 1.3 2004/10/01 16:33:39 ivos
! cosmetics.
!
! Revision 1.2 2004/10/01 16:09:11 ivos
! Replaced REAL(ppm_kind_double) :: t0 with REAL(MK) t0.
!
! Revision 1.1 2004/09/23 09:56:40 kotsalie
! Mg new version
!
! Revision 1.4 2004/07/26 13:49:19 ivos
! Removed Routines sections from the header comment.
!
! Revision 1.3 2004/07/26 11:59:40 ivos
! Fixes to make it compile.
!
! Revision 1.2 2004/07/26 07:42:39 ivos
! Changed to single-interface modules and adapted all USE statements.
!
! Revision 1.1 2004/06/29 14:36:49 kotsalie
! Commiting multigrid for further use
!
!-------------------------------------------------------------------------
! Copyright (c) 2012 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
! Center for Fluid Dynamics (DTU)
!
!
! This file is part of the Parallel Particle Mesh Library (PPM).
!
! PPM is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License
! as published by the Free Software Foundation, either
! version 3 of the License, or (at your option) any later
! version.
!
! PPM is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! and the GNU Lesser General Public License along with PPM. If not,
! see <http://www.gnu.org/licenses/>.
!
! Parallel Particle Mesh Library (PPM)
! ETH Zurich
! CH-8092 Zurich, Switzerland
!-------------------------------------------------------------------------
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_finalize_2d_sca_s(info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_finalize_2d_sca_d(info)
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_finalize_3d_sca_s(info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_finalize_3d_sca_d(info)
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_finalize_2d_vec_s(info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_finalize_2d_vec_d(info)
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
SUBROUTINE ppm_mg_finalize_3d_vec_s(info)
#elif __KIND == __DOUBLE_PRECISION
SUBROUTINE ppm_mg_finalize_3d_vec_d(info)
#endif
#endif
#endif
!-------------------------------------------------------------------------
! Includes
!-------------------------------------------------------------------------
#include "ppm_define.h"
!-------------------------------------------------------------------------
! Modules
!-------------------------------------------------------------------------
USE ppm_module_data
USE ppm_module_data_mg
USE ppm_module_mg_alloc
USE ppm_module_substart
USE ppm_module_substop
USE ppm_module_error
USE ppm_module_alloc
USE ppm_module_write
IMPLICIT NONE
#if __KIND == __SINGLE_PRECISION | __KIND == __SINGLE_PRECISION_COMPLEX
INTEGER, PARAMETER :: MK = ppm_kind_single
#else
INTEGER, PARAMETER :: MK = ppm_kind_double
#endif
!-------------------------------------------------------------------------
! Arguments
!-------------------------------------------------------------------------
INTEGER, INTENT(INOUT) :: info
!-------------------------------------------------------------------------
! Local variables
!-------------------------------------------------------------------------
INTEGER, DIMENSION(1) :: lda1
INTEGER, DIMENSION(2) :: lda2
INTEGER, DIMENSION(3) :: lda3
INTEGER, DIMENSION(4) :: lda4
INTEGER, DIMENSION(5) :: lda5
INTEGER :: iopt,i
INTEGER :: istat,j
REAL(MK) :: t0
CHARACTER(LEN=ppm_char) :: mesg
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
TYPE(mg_field_2d_sca_s),DIMENSION(:,:),POINTER :: mgfield
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_2d_sca_d),DIMENSION(:,:),POINTER :: mgfield
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
TYPE(mg_field_3d_sca_s),DIMENSION(:,:),POINTER :: mgfield
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_3d_sca_d),DIMENSION(:,:),POINTER :: mgfield
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
TYPE(mg_field_2d_vec_s),DIMENSION(:,:),POINTER :: mgfield
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_2d_vec_d),DIMENSION(:,:),POINTER :: mgfield
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
TYPE(mg_field_3d_vec_s),DIMENSION(:,:),POINTER :: mgfield
#elif __KIND == __DOUBLE_PRECISION
TYPE(mg_field_3d_vec_d),DIMENSION(:,:),POINTER :: mgfield
#endif
#endif
#endif
!-------------------------------------------------------------------------
! Externals
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Initialise
!-------------------------------------------------------------------------
CALL substart('ppm_mg_finalize',t0,info)
lda1(1)=0
lda2(1:2) = 0
lda3(1:3) = 0
lda4(1:4) = 0
lda5(1:5) = 0
!-------------------------------------------------------------------------
!Definition of necessary variables and allocation of arrays
!-------------------------------------------------------------------------
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
mgfield=>mgfield_2d_sca_s
#elif __KIND == __DOUBLE_PRECISION
mgfield=>mgfield_2d_sca_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
mgfield=>mgfield_3d_sca_s
#elif __KIND == __DOUBLE_PRECISION
mgfield=>mgfield_3d_sca_d
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
mgfield=>mgfield_2d_vec_s
#elif __KIND == __DOUBLE_PRECISION
mgfield=>mgfield_2d_vec_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
mgfield=>mgfield_3d_vec_s
#elif __KIND == __DOUBLE_PRECISION
mgfield=>mgfield_3d_vec_d
#endif
#endif
#endif
!-------------------------------------------------------------------------
! Deallocate global arrays (from ppm_module_multigrid)
!-------------------------------------------------------------------------
istat = 0
iopt = ppm_param_dealloc
CALL ppm_alloc(start,lda3,iopt,info)
istat=istat+info
CALL ppm_alloc(lboundary,lda2,iopt,info)
istat=istat+info
CALL ppm_alloc(max_node,lda2,iopt,info)
istat=istat+info
#if __DIM == __SFIELD
CALL ppm_alloc(bcdef_sca,lda1,iopt,info)
istat=istat+info
#elif __DIM == __VFIELD
CALL ppm_alloc(bcdef_vec,lda1,iopt,info)
istat=istat+info
#endif
CALL ppm_alloc(ghostsize,lda1,iopt,info)
istat=istat+info
CALL ppm_alloc(factor,lda1,iopt,info)
istat=istat+info
CALL ppm_alloc(mg_meshid,lda1,iopt,info)
istat=istat+info
CALL ppm_mg_alloc(mgfield,lda2,iopt,info)
istat = istat +info
IF (istat .NE. 0) THEN
WRITE(mesg,'(A,I3,A)') 'for ',istat,' mgr arrays.Pble memory leak.'
info = ppm_error_error
CALL ppm_error(ppm_err_dealloc,'ppm_mg_finalize',mesg,__LINE__,&
& info)
GOTO 9999
ENDIF
!-------------------------------------------------------------------------
! Return
!-------------------------------------------------------------------------
9999 CONTINUE
CALL substop('ppm_mg_finalize',t0,info)
RETURN
#if __DIM == __SFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_finalize_2d_sca_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_finalize_2d_sca_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_finalize_3d_sca_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_finalize_3d_sca_d
#endif
#endif
#elif __DIM == __VFIELD
#if __MESH_DIM == __2D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_finalize_2d_vec_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_finalize_2d_vec_d
#endif
#elif __MESH_DIM == __3D
#if __KIND == __SINGLE_PRECISION
END SUBROUTINE ppm_mg_finalize_3d_vec_s
#elif __KIND == __DOUBLE_PRECISION
END SUBROUTINE ppm_mg_finalize_3d_vec_d
#endif
#endif
#endif
src/mg/ppm_mg_init.f deleted 100644 → 0
+ 0
− 1158
View file @ c54e8f1f
This diff is collapsed.
src/mg/ppm_mg_prolong.f deleted 100644 → 0
+ 0
− 564
View file @ c54e8f1f
This diff is collapsed.
src/mg/ppm_mg_res_coarse.f deleted 100644 → 0
+ 0
− 481
View file @ c54e8f1f
This diff is collapsed.
src/mg/ppm_mg_res_fine.f deleted 100644 → 0
+ 0
− 465
View file @ c54e8f1f
This diff is collapsed.
src/mg/ppm_mg_restrict.f deleted 100644 → 0
+ 0
− 775
View file @ c54e8f1f
This diff is collapsed.
src/mg/ppm_mg_smooth_coarse.f deleted 100644 → 0
+ 0
− 1112
View file @ c54e8f1f
This diff is collapsed.
src/mg/ppm_mg_smooth_fine.f deleted 100644 → 0
+ 0
− 887
View file @ c54e8f1f
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment