From 6b6f3c7728419d520feb85a0c3d9be7ea6382d8c Mon Sep 17 00:00:00 2001
From: jtra <jtra@7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9>
Date: Tue, 21 Jun 2011 11:31:14 +0000
Subject: [PATCH] included freespace solution to the poisson equation included
spectral curl - but ONLY for the periodic case cleanup, maybe small bugs
git-svn-id: https://ppm.inf.ethz.ch/svn/ppmnumerics/branches/ngtopo/libppmnumerics@1024 7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9
---
...isson_init_predef.f => ppm_poisson_init.f} | 110 ++----
src/poisson/ppm_poisson_solve.f | 363 ++++--------------
src/ppm_module_poisson.f | 6 +-
3 files changed, 111 insertions(+), 368 deletions(-)
rename src/poisson/{ppm_poisson_init_predef.f => ppm_poisson_init.f} (85%)
diff --git a/src/poisson/ppm_poisson_init_predef.f b/src/poisson/ppm_poisson_init.f
similarity index 85%
rename from src/poisson/ppm_poisson_init_predef.f
rename to src/poisson/ppm_poisson_init.f
index 77e2d98..38caaaa 100644
--- a/src/poisson/ppm_poisson_init_predef.f
+++ b/src/poisson/ppm_poisson_init.f
@@ -1,5 +1,5 @@
!-------------------------------------------------------------------------
- ! Subroutine : ppm_poisson_init_predef.f90
+ ! Subroutine : ppm_poisson_init.f90
!-------------------------------------------------------------------------
! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
! Center for Fluid Dynamics (DTU)
@@ -50,7 +50,7 @@
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_pois_fre - Poisson equation, freespace boundaries
!!!ppm_poisson_grn_reprojec - Do vorticity reprojection to kill divergence
!!!
!!!Eventually this should also accept custom Greens function
@@ -64,7 +64,7 @@
!!!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_sp (not fully implemented)
+ !!! * ppm_poisson_drv_curl_sp (only for periodic BC)
!!! * ppm_poisson_drv_grad_sp (not implemented)
!!! * ppm_poisson_drv_lapl_sp (not implemented)
!!! * ppm_poisson_drv_div_sp (not implemented)
@@ -111,7 +111,7 @@
!-------------------------------------------------------------------------
! Initialise routine
!-------------------------------------------------------------------------
- CALL substart('ppm_poisson_init_predef',t0,info)
+ CALL substart('ppm_poisson_init',t0,info)
!-------------------------------------------------------------------------
! Investigate optional arguments, setup routine accordingly
@@ -130,7 +130,7 @@
!-------------------------------------------------------------------------
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)
+ CALL ppm_write(ppm_rank,'ppm_poisson_init','Failed to get topology.',isub)
GOTO 9999
ENDIF
!nsubs = topology%nsublist
@@ -173,10 +173,6 @@
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
!-------------------------------------------------------------------------
@@ -192,33 +188,29 @@
& 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
+ IF (ppmpoisson%case .EQ. ppm_poisson_grn_pois_per) THEN
+ ppmpoisson%normkx = &
+ & 2.0_MK*PI/(topology%max_physd(1)-topology%min_physd(1))
+ ppmpoisson%normky = &
+ & 2.0_MK*PI/(topology%max_physd(2)-topology%min_physd(2))
+ ppmpoisson%normkz = &
+ & 2.0_MK*PI/(topology%max_physd(3)-topology%min_physd(3))
+ ELSE IF (ppmpoisson%case .EQ. ppm_poisson_grn_pois_fre) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_init', &
+ & 'WARNING: Spectral curl is not fully implemented in Freespace.',isub)
+ ppmpoisson%normkx = &
+ & 2.0_MK*PI/((topology%max_physd(1)-topology%min_physd(1))*2.0_MK)
+ ppmpoisson%normky = &
+ & 2.0_MK*PI/((topology%max_physd(2)-topology%min_physd(2))*2.0_MK)
+ ppmpoisson%normkz = &
+ & 2.0_MK*PI/((topology%max_physd(3)-topology%min_physd(3))*2.0_MK)
+ ENDIF
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)
+ CALL ppm_write(ppm_rank,'ppm_poisson_init','Undefined derivation input.',isub)
GOTO 9999
ENDIF
ELSE
@@ -228,7 +220,7 @@
!-------------------------------------------------------------------------
- ! Create new slab topology !@sofar periodic
+ ! Create new slab topology
!-------------------------------------------------------------------------
ttopoid = 0
tmeshid = -1
@@ -241,8 +233,6 @@
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,&
@@ -251,7 +241,7 @@
& __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)
+ CALL ppm_write(ppm_rank,'ppm_poisson_init','Failed to create xy-topology.',isub)
GOTO 9999
ENDIF
ppmpoisson%topoidxy = ttopoid
@@ -266,16 +256,14 @@
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)
+ CALL ppm_write(ppm_rank,'ppm_poisson_init','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
+ ! Create new pencil topology
!-------------------------------------------------------------------------
ttopoid = 0
tmeshid = -1
@@ -286,17 +274,14 @@
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)
+ CALL ppm_write(ppm_rank,'ppm_poisson_init','Failed to create z-topology.',isub)
GOTO 9999
ENDIF
ppmpoisson%topoidz = ttopoid
@@ -307,19 +292,19 @@
!-------------------------------------------------------------------------
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)
+ CALL ppm_write(ppm_rank,'ppm_poisson_init','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)
+ !!CALL ppm_write(ppm_rank,'ppm_poisson_init','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)
+ CALL ppm_write(ppm_rank,'ppm_poisson_init','Failed to get z topology.',isub)
GOTO 9999
ENDIF
@@ -389,12 +374,8 @@
ENDDO
!-------------------------------------------------------------------------
- ! Allocate real and complex xy slabs
+ ! Allocate 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)
@@ -434,11 +415,6 @@
& 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)
@@ -474,10 +450,13 @@
! 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?
+ ! Scaling the spectral coefficients...
+ ! one minus due to (i*k)^2 and another due to the Poisson equation
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))
+ & 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)
@@ -524,18 +503,11 @@
! First initialise the real Greens function
!@alternatively this could come from as input
!-----------------------------------------------------------------------
- !@write(*,*) 'what the fuck?'
- !there should NOT be a minus here since this Greens function takes
+ !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 greens 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)
@@ -638,13 +610,11 @@
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
+ ppmpoisson%fldgrnc(i,j,k,isub) = ppmpoisson%fldzc2(1,i,j,k,isub)
ENDDO
ENDDO
ENDDO
ENDDO
-
- !@Deallocate field green r NOPE!
END IF
!-------------------------------------------------------------------------
@@ -659,7 +629,7 @@
! Return
!-------------------------------------------------------------------------
9999 CONTINUE
- CALL substop('ppm_poisson_init_predef',t0,info)
+ CALL substop('ppm_poisson_init',t0,info)
RETURN
END SUBROUTINE __ROUTINE
diff --git a/src/poisson/ppm_poisson_solve.f b/src/poisson/ppm_poisson_solve.f
index e9e0382..7488423 100644
--- a/src/poisson/ppm_poisson_solve.f
+++ b/src/poisson/ppm_poisson_solve.f
@@ -54,16 +54,13 @@
INTEGER :: i,j,k
INTEGER :: info2
INTEGER :: presentcase
- REAL(__PREC) :: wdotk
+ REAL(__PREC) :: wdotgreenr
+ COMPLEX(__PREC) :: wdotgreenc
INTEGER :: gi,gj,gk
- REAL(__PREC) :: kx,ky,kz
- REAL(__PREC) :: phix,phiy,phiz
+ COMPLEX(__PREC) :: kx,ky,kz
+ COMPLEX(__PREC) :: phix,phiy,phiz
REAL(__PREC) :: normfac
-#ifndef __NOPE
-INTEGER :: trank !@
-trank =0
-#endif
!-------------------------------------------------------------------------
! Initialise routine
!-------------------------------------------------------------------------
@@ -144,94 +141,6 @@ trank =0
& ppmpoisson%fldxyr, ppmpoisson%fldxyc, &
& info)
-#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
!-----------------------------------------------------------------------
! Map to the pencils (Z)
@@ -283,67 +192,6 @@ trank =0
& ppmpoisson%fldzc1, ppmpoisson%fldzc2, &
& info)
-#ifdef __NOPE
- if (ppm_rank .EQ. trank) THEN
- write(*,*)
- write(*,*) !ppmpoisson%fldzc1 = 0.0_MK
- write(*,*) 'fftx123', ppm_rank
- DO isub=1,ppmpoisson%nsublistz
- isubl=ppmpoisson%isublistz(isub)
- DO k=1,ppmpoisson%ndataz(3,isubl)-1
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataz(1,isubl)-1
- DO j=1,ppmpoisson%ndataz(2,isubl)-1
- write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldzc2(1,i,j,k,isub),')'
- END DO
- write(*,*)
- END DO
- END DO
- END DO
- write(*,*) 'ffty123', ppm_rank
- DO isub=1,ppmpoisson%nsublistz
- isubl=ppmpoisson%isublistz(isub)
- DO k=1,ppmpoisson%ndataz(3,isubl)-1
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataz(1,isubl)-1
- DO j=1,ppmpoisson%ndataz(2,isubl)-1
- write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldzc2(2,i,j,k,isub),')'
- END DO
- write(*,*)
- END DO
- END DO
- END DO
- write(*,*) 'fftz123', ppm_rank
- DO isub=1,ppmpoisson%nsublistz
- isubl=ppmpoisson%isublistz(isub)
- DO k=1,ppmpoisson%ndataz(3,isubl)-1
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataz(1,isubl)-1
- DO j=1,ppmpoisson%ndataz(2,isubl)-1
- write(*,'(A,E12.4,E12.4,A,$)') ' (',ppmpoisson%fldzc2(3,i,j,k,isub),')'
- END DO
- write(*,*)
- END DO
- END DO
- END DO
-
- write(*,*)
- write(*,*)
- write(*,*) 'green', ppm_rank
- DO isub=1,ppmpoisson%nsublistz
- isubl=ppmpoisson%isublistz(isub)
- DO k=1,ppmpoisson%ndataz(3,isubl)-1
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataz(1,isubl)-1
- DO j=1,ppmpoisson%ndataz(2,isubl)-1
- write(*,'(E12.4,$)') REAL(ppmpoisson%fldgrnr(i,j,k,isub))
- END DO
- write(*,*)
- END DO
- END DO
- END DO
- ENDIF
-#endif
!-----------------------------------------------------------------------
! Apply the periodic Greens function
@@ -384,90 +232,109 @@ trank =0
ENDDO
ENDDO
!-----------------------------------------------------------------------
- ! Spectral derivatives
- ! normkx, etc contains 2pi/Lx
+ ! Vorticity re-projection
!-----------------------------------------------------------------------
- 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!'
+ 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)
+ !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)
+ 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 = CMPLX(0.0_MK,gj*ppmpoisson%normky,MK)
+ 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 = CMPLX(0.0_MK,gi*ppmpoisson%normkx,MK)
+ kx = REAL(gi,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)
+ !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
+
+ IF (ppmpoisson%case .EQ. ppm_poisson_grn_pois_fre) THEN
+ wdotgreenc = (ppmpoisson%fldzc2(1,i,j,k,isub) * kx + &
+ & ppmpoisson%fldzc2(2,i,j,k,isub) * ky + &
+ & ppmpoisson%fldzc2(3,i,j,k,isub) * kz) * &
+ & ppmpoisson%fldgrnc( i,j,k,isub)
+ ppmpoisson%fldzc2(1,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(1,i,j,k,isub) - wdotgreenc*kx)
+ ppmpoisson%fldzc2(2,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(2,i,j,k,isub) - wdotgreenc*ky)
+ ppmpoisson%fldzc2(3,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(3,i,j,k,isub) - wdotgreenc*kz)
+ ELSE IF (ppmpoisson%case .EQ. ppm_poisson_grn_pois_per) THEN
+ wdotgreenr = (ppmpoisson%fldzc2(1,i,j,k,isub) * kx + &
+ & ppmpoisson%fldzc2(2,i,j,k,isub) * ky + &
+ & ppmpoisson%fldzc2(3,i,j,k,isub) * kz) * &
+ & ppmpoisson%fldgrnr( i,j,k,isub)
+ ppmpoisson%fldzc2(1,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(1,i,j,k,isub) - wdotgreenr*kx)
+ ppmpoisson%fldzc2(2,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(2,i,j,k,isub) - wdotgreenr*ky)
+ ppmpoisson%fldzc2(3,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(3,i,j,k,isub) - wdotgreenr*kz)
+ ENDIF
- !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
+ ! Spectral derivatives
+ ! normkx, etc contains 2pi/Lx
!-----------------------------------------------------------------------
- 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)
+ 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
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)
+ 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 = REAL(gj,MK)
+ 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 = REAL(gi,MK)
+ kx = CMPLX(0.0_MK,gi*ppmpoisson%normkx,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
+ phix = ppmpoisson%fldzc2(1,i,j,k,isub)
+ phiy = ppmpoisson%fldzc2(2,i,j,k,isub)
+ phiz = ppmpoisson%fldzc2(3,i,j,k,isub)
- 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
+ ppmpoisson%fldzc2(1,i,j,k,isub) = (ky*phiz-kz*phiy)
+ ppmpoisson%fldzc2(2,i,j,k,isub) = (kz*phix-kx*phiz)
+ ppmpoisson%fldzc2(3,i,j,k,isub) = (kx*phiy-ky*phix)
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
+
!-----------------------------------------------------------------------
! IFFT pencil (Z)
!-----------------------------------------------------------------------
@@ -476,51 +343,6 @@ trank =0
& ppmpoisson%fldzc2, ppmpoisson%fldzc1, &
& info)
-#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)
@@ -570,51 +392,6 @@ trank =0
& ppmpoisson%fldxyc, ppmpoisson%fldxyr, &
& info)
-#ifdef __NOPE
- if (ppm_rank .EQ. trank) THEN
- write(*,*)
- write(*,*)
- write(*,*) 'ifftx', ppm_rank
- DO isub=1,ppmpoisson%nsublistxy
- isubl=ppmpoisson%isublistxy(isub)
- DO k=1,ppmpoisson%ndataxy(3,isubl)-1
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataxy(1,isubl)-1
- DO j=1,ppmpoisson%ndataxy(2,isubl)-1
- write(*,'(E12.4,$)') ppmpoisson%fldxyr(1,i,j,k,isub)
- END DO
- write(*,*)
- END DO
- END DO
- END DO
- write(*,*) 'iffty', ppm_rank
- DO isub=1,ppmpoisson%nsublistxy
- isubl=ppmpoisson%isublistxy(isub)
- DO k=1,ppmpoisson%ndataxy(3,isubl)-1
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataxy(1,isubl)-1
- DO j=1,ppmpoisson%ndataxy(2,isubl)-1
- write(*,'(E12.4,$)') ppmpoisson%fldxyr(2,i,j,k,isub)
- END DO
- write(*,*)
- END DO
- END DO
- END DO
- write(*,*) 'ifftz', ppm_rank
- DO isub=1,ppmpoisson%nsublistxy
- isubl=ppmpoisson%isublistxy(isub)
- DO k=1,ppmpoisson%ndataxy(3,isubl)-1
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataxy(1,isubl)-1
- DO j=1,ppmpoisson%ndataxy(2,isubl)-1
- write(*,'(E12.4,$)') ppmpoisson%fldxyr(3,i,j,k,isub)
- END DO
- write(*,*)
- END DO
- END DO
- END DO
- ENDIF
-#endif
!-----------------------------------------------------------------------
! Map back to standard topology (XYZ)
@@ -718,15 +495,11 @@ trank =0
! 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
+ IF (ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_fd2) 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
+ IF (ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_fd4) THEN
CALL ppm_poisson_fd(topoid,meshid,ppmpoisson%drv_vr,fieldout,&
& ppm_poisson_drv_curl_fd4,info)
ENDIF
diff --git a/src/ppm_module_poisson.f b/src/ppm_module_poisson.f
index fd0408d..e044d1a 100644
--- a/src/ppm_module_poisson.f
+++ b/src/ppm_module_poisson.f
@@ -176,7 +176,7 @@
INTEGER,PARAMETER :: ppm_poisson_grn_reprojec =3
INTERFACE ppm_poisson_init
- MODULE PROCEDURE ppm_poisson_init_predef
+ MODULE PROCEDURE ppm_poisson_init
END INTERFACE
INTERFACE ppm_poisson_solve
@@ -199,8 +199,8 @@
#define __DIM 3
#define __ZEROSI (/0,0,0/)
#define __NCOM 3
-#define __ROUTINE ppm_poisson_init_predef
-#include "poisson/ppm_poisson_init_predef.f"
+#define __ROUTINE ppm_poisson_init
+#include "poisson/ppm_poisson_init.f"
#undef __ROUTINE
#define __ROUTINE ppm_poisson_solve
#include "poisson/ppm_poisson_solve.f"
--
GitLab