From caf17f98faa663d0e48f0bfe0b2a7cfa798dc2da Mon Sep 17 00:00:00 2001
From: jtra <jtra@7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9>
Date: Sun, 26 Jun 2011 08:32:44 +0000
Subject: [PATCH] new features: spectral derivatives(curl) for periodic greens
function vorticity reprojection for periodic domains freespace solution to
poisson equation
git-svn-id: https://ppm.inf.ethz.ch/svn/ppmnumerics/branches/ngtopo/libppmnumerics@1025 7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9
---
src/poisson/ppm_poisson_init.f | 35 +++++++++------
src/poisson/ppm_poisson_solve.f | 77 ++++++++++++++++++---------------
2 files changed, 63 insertions(+), 49 deletions(-)
diff --git a/src/poisson/ppm_poisson_init.f b/src/poisson/ppm_poisson_init.f
index 38caaaa..b48b9ab 100644
--- a/src/poisson/ppm_poisson_init.f
+++ b/src/poisson/ppm_poisson_init.f
@@ -190,22 +190,9 @@
& LBOUND(fieldout,5):UBOUND(fieldout,5)))
ELSE IF (derive .EQ. ppm_poisson_drv_curl_sp) THEN
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
+ 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
@@ -216,6 +203,26 @@
ELSE
ppmpoisson%derivatives = ppm_poisson_drv_none
ENDIF
+ !-------------------------------------------------------------------------
+ ! Create spectral scaling components always. Just in case some
+ ! reprojection comes up
+ ! The conditionals need to be for not just the Poisson equation
+ !-------------------------------------------------------------------------
+ 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
+ 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
diff --git a/src/poisson/ppm_poisson_solve.f b/src/poisson/ppm_poisson_solve.f
index 7488423..12ab83c 100644
--- a/src/poisson/ppm_poisson_solve.f
+++ b/src/poisson/ppm_poisson_solve.f
@@ -54,8 +54,7 @@
INTEGER :: i,j,k
INTEGER :: info2
INTEGER :: presentcase
- REAL(__PREC) :: wdotgreenr
- COMPLEX(__PREC) :: wdotgreenc
+ COMPLEX(__PREC) :: helmholtzpic
INTEGER :: gi,gj,gk
COMPLEX(__PREC) :: kx,ky,kz
COMPLEX(__PREC) :: phix,phiy,phiz
@@ -236,23 +235,29 @@
!-----------------------------------------------------------------------
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%case .EQ. ppm_poisson_grn_pois_fre) THEN
+ normfac = 1.0_MK/ REAL((ppmpoisson%nmz(1))* & !vertex
+ & (ppmpoisson%nmz(2))* &
+ & (ppmpoisson%nmz(3)),MK)
+ ELSE IF (ppmpoisson%case .EQ. ppm_poisson_grn_pois_per) THEN
+ normfac = 1.0_MK/ REAL((ppmpoisson%nmz(1)-1)* & !vertex
+ & (ppmpoisson%nmz(2)-1)* &
+ & (ppmpoisson%nmz(3)-1),MK)
+ ENDIF
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 = REAL(gk,MK)*ppmpoisson%normkz
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 = REAL(gj,MK)*ppmpoisson%normky
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 = REAL(gi,MK)*ppmpoisson%normkx
!IF (gi .EQ. 0 .AND. gj .EQ. 0 .AND. gk .EQ. 0) THEN
!wdotk = 0.0_mk
@@ -271,27 +276,27 @@
!& (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)
+ helmholtzpic = (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) = & !@formatting
+ & (ppmpoisson%fldzc2(1,i,j,k,isub)*normfac - helmholtzpic*kx)
+ ppmpoisson%fldzc2(2,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(2,i,j,k,isub)*normfac - helmholtzpic*ky)
+ ppmpoisson%fldzc2(3,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(3,i,j,k,isub)*normfac - helmholtzpic*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)
+ helmholtzpic = (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) = & !@formatting
+ & (ppmpoisson%fldzc2(1,i,j,k,isub)*normfac - helmholtzpic*kx)
+ ppmpoisson%fldzc2(2,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(2,i,j,k,isub)*normfac - helmholtzpic*ky)
+ ppmpoisson%fldzc2(3,i,j,k,isub) = &
+ & (ppmpoisson%fldzc2(3,i,j,k,isub)*normfac - helmholtzpic*kz)
ENDIF
ENDDO
@@ -495,13 +500,15 @@
! 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) 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) THEN
- CALL ppm_poisson_fd(topoid,meshid,ppmpoisson%drv_vr,fieldout,&
- & ppm_poisson_drv_curl_fd4,info)
+ IF (presentcase .NE. ppm_poisson_grn_reprojec) 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) THEN
+ CALL ppm_poisson_fd(topoid,meshid,ppmpoisson%drv_vr,fieldout,&
+ & ppm_poisson_drv_curl_fd4,info)
+ ENDIF
ENDIF
!-------------------------------------------------------------------------
--
GitLab