From 83ac556d0e30092566e6fe0e0895d7959a495a79 Mon Sep 17 00:00:00 2001
From: jtra <jtra@7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9>
Date: Wed, 15 Jun 2011 14:35:05 +0000
Subject: [PATCH] included freespace solution to poisson solver
git-svn-id: https://ppm.inf.ethz.ch/svn/ppmnumerics/branches/ngtopo/libppmnumerics@888 7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9
---
src/poisson/ppm_poisson_extrapolateghost.f | 32 ++-
src/poisson/ppm_poisson_fd.f | 43 +++-
src/poisson/ppm_poisson_init_predef.f | 65 ++++--
src/poisson/ppm_poisson_solve.f | 249 ++++++++-------------
src/ppm_define.h | 1 +
src/ppm_module_poisson.f | 4 +
6 files changed, 203 insertions(+), 191 deletions(-)
diff --git a/src/poisson/ppm_poisson_extrapolateghost.f b/src/poisson/ppm_poisson_extrapolateghost.f
index c2cc564..980e624 100644
--- a/src/poisson/ppm_poisson_extrapolateghost.f
+++ b/src/poisson/ppm_poisson_extrapolateghost.f
@@ -1,11 +1,23 @@
!-------------------------------------------------------------------------
- ! ppm_poisson_extrapolateghost.f90
+ ! Subroutine : ppm_poisson_extrapolateghost.f90
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
!-------------------------------------------------------------------------
-!!#define __ROUTINE ppm_poisson_extrapolateghost_vr
-!!#define __DIM 3
-!!#define __NCOM 3
-!!#define __ZEROSI (/0,0,0/)
SUBROUTINE __ROUTINE(topoid,meshid,field,nextra,nbase,gstw,info)
+ !!! This routine extrapolates field values of field with topology and mesh
+ !!! id topoid,meshid respectively into the ghost layer of width gstw.
+ !!! nbase points are used to extrapolate into nextra points.
+ !!!
+ !!! [NOTE]
+ !!! Presently extrapolation can only be done to 1 or two points into the
+ !!! ghostlayer always based on 4 points (fourth order spatial convergence)
+ !!! A general nbase,nextra extrapolation can be implemented vi solution
+ !!! of a small linear system of equations. This has not been done.
+ !!! This routine is in need of loop unrollling in particular for
+ !!! typical choices of nbase,nextra pairs. Extrapolation is necessary for
+ !!! e.g. freespace FD curl of stream function
USE ppm_module_topo_get
@@ -51,19 +63,15 @@
!-------------------------------------------------------------------------
! Compare the number of points to extrapolate to the ghost layer width
!-------------------------------------------------------------------------
- IF (iextra .GT. gstw(1) .OR. &
- & iextra .GT. gstw(2) .OR. &
- & iextra .GT. gstw(3)) THEN
+ IF (nextra .GT. gstw(1) .OR. &
+ & nextra .GT. gstw(2) .OR. &
+ & nextra .GT. gstw(3)) THEN
CALL ppm_write(ppm_rank,'ppm_poisson_extrapolateghost',&
& 'The points to extrapolate exceeds the ghost layer.',info)
info = -1
GOTO 9999
ENDIF
- !-------------------------------------------------------------------------
- !@ Check what has been implemented in this routine so far
- !-------------------------------------------------------------------------
-
!-------------------------------------------------------------------------
! Determine weights
!-------------------------------------------------------------------------
diff --git a/src/poisson/ppm_poisson_fd.f b/src/poisson/ppm_poisson_fd.f
index ebc7415..383f9f4 100644
--- a/src/poisson/ppm_poisson_fd.f
+++ b/src/poisson/ppm_poisson_fd.f
@@ -1,13 +1,28 @@
!-------------------------------------------------------------------------
- ! ppm_poisson_fd.f90
+ ! Subroutine : ppm_poisson_fd.f90
!-------------------------------------------------------------------------
- !@ TODO: Somewhere check if fieldin is equal to fieldout and give a warning
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
!-------------------------------------------------------------------------
- !!#define __ROUTINE ppm_poisson_fd
- !!#define __DIM 3
- !!#define __NCOM 3
- !!#define __ZEROSI (/0,0,0/)
SUBROUTINE __ROUTINE(topoid,meshid,fieldin,fieldout,dtype,info)
+ !!! This routine computes the finite difference gradients, curl etc of
+ !!! fieldin and outputs to fieldout. Both in and out fields are on
+ !!! the mesh with meshid belonging to the topology with id topoid.
+ !!! The finite difference to be carried out is determined by dtype which
+ !!! must be one of the following:
+ !!! * ppm_poisson_drv_curl_fd2
+ !!! * ppm_poisson_drv_grad_fd2 (not implemented yet)
+ !!! * ppm_poisson_drv_lapl_fd2 (not implemented yet)
+ !!! * ppm_poisson_drv_div_fd2 (not implemented yet)
+ !!! * ppm_poisson_drv_curl_fd4
+ !!! * ppm_poisson_drv_grad_fd4 (not implemented yet)
+ !!! * ppm_poisson_drv_lapl_fd4 (not implemented yet)
+ !!! * ppm_poisson_drv_div_fd4 (not implemented yet)
+ !!!
+ !!! [NOTE] fieldin and fieldout must NOT be the same array. A check
+ !!! should be added.
+ !@ TODO: Somewhere check if fieldin is equal to fieldout and give a warning
USE ppm_module_topo_get
@@ -16,12 +31,25 @@
! Arguments
!-------------------------------------------------------------------------
INTEGER, INTENT(IN) :: topoid
+ !!! ID of the topology
INTEGER, INTENT(IN) :: meshid
+ !!! Mesh ID
REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldin
+ !!! Input field data
REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldout
- !INTEGER,DIMENSION(__DIM),INTENT(IN) :: gstw
+ !!! Output field data
INTEGER, INTENT(IN) :: dtype
+ !!! Derivation type. Can be one of the types:
+ !!! * ppm_poisson_drv_curl_fd2
+ !!! * ppm_poisson_drv_grad_fd2 (not implemented yet)
+ !!! * ppm_poisson_drv_lapl_fd2 (not implemented yet)
+ !!! * ppm_poisson_drv_div_fd2 (not implemented yet)
+ !!! * ppm_poisson_drv_curl_fd4
+ !!! * ppm_poisson_drv_grad_fd4 (not implemented yet)
+ !!! * ppm_poisson_drv_lapl_fd4 (not implemented yet)
+ !!! * ppm_poisson_drv_div_fd4 (not implemented yet)
INTEGER, INTENT(OUT) :: info
+ !!! Return status, 0 upon succes
!-------------------------------------------------------------------------
! Local variables
@@ -40,7 +68,6 @@
!-------------------------------------------------------------------------
CALL substart('ppm_poisson_fd',t0,info)
- !@ perhaps ensure that fieldin .NE. fieldout
!-------------------------------------------------------------------------
! Get topology and mesh values
!-------------------------------------------------------------------------
diff --git a/src/poisson/ppm_poisson_init_predef.f b/src/poisson/ppm_poisson_init_predef.f
index 5baa4c8..104afef 100644
--- a/src/poisson/ppm_poisson_init_predef.f
+++ b/src/poisson/ppm_poisson_init_predef.f
@@ -1,16 +1,29 @@
!-------------------------------------------------------------------------
- ! ppm_poisson_init_predef.f90
+ ! Subroutine : ppm_poisson_init_predef.f90
!-------------------------------------------------------------------------
- ! Routine to initialise the convolution using build in Greens function
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
!-------------------------------------------------------------------------
-!!#define __PREC ppm_kind_double
-!!#define __CMPLXDEF DCMPLX
-!!#define __DIM 3
-!!#define __NCOM 3
-!!#define __ZEROSI (/0,0,0/)
-!!#define __ROUTINE ppm_poisson_init_predef
SUBROUTINE __ROUTINE(topoid,meshid,ppmpoisson,fieldin,fieldout,green,info&
&,bc,derive)
+ !!! Routine to initialise Green's function solution of the Poisson
+ !!! equation. green is the flag defining which Green's function to use:
+ !!! * ppm_poisson_grn_pois_per - Poisson equation, periodic boundaries
+ !!! * ppm_poisson_grn_pois_fre - Poisson equation, freespace boundaries (not implemented)
+ !!! * ppm_poisson_grn_reprojec - Do vorticity reprojection to kill divergence
+ !!! Eventually the routine should be overloaded to accept custom Green's
+ !!! functions such that more general convolutions can be performed.
+ !!! green should be expanded to include more buildin Green's functions.
+ !!!
+ !!! The routine should accept an optional flag to toggle deallocation of
+ !!! work arrays between calls to ppm_poisson_solve
+ !!!
+ !!! [NOTE]
+ !!! When meshid > 1 works with the mapping the memory consumption can be
+ !!! halved. Sketches have been implemented.
+ !!! This routine should also be renamed to _init
+
USE ppm_module_mktopo
USE ppm_module_topo_get
@@ -24,33 +37,49 @@
! Arguments
!-------------------------------------------------------------------------
INTEGER, INTENT(IN) :: topoid
+ !!! Topology ID
INTEGER, INTENT(IN) :: meshid
+ !!! Mesh ID
TYPE(ppm_poisson_plan),INTENT(INOUT) :: ppmpoisson
- !@ strictly speaking fieldin is not being used in the init routine
+ !!! The PPM Poisson plan type (inspired by the FFTW plan)
REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldin
+ !!! Input data field
+ !@ strictly speaking fieldin is not being used in the init routine
REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldout
- !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: green
+ !!! Output data field
INTEGER, INTENT(IN) :: green
!!!flag to select build-in Greens functions:
!!!ppm_poisson_grn_pois_per - Poisson equation, periodic boundaries
!!!ppm_poisson_grn_pois_fre - Poisson equation, freespace boundaries (not implemented)
!!!ppm_poisson_grn_reprojec - Do vorticity reprojection to kill divergence
+ !!!
+ !!!Eventually this should also accept custom Green's function
INTEGER, INTENT(OUT) :: info
INTEGER,INTENT(IN),OPTIONAL :: bc
!!!boundary condition for the convolution. Can be on of the following:
!!!ppm_poisson_grn_pois_per, ppm_poisson_grn_pois_fre.
!!!One could argue that this is redundant in the build-in case
- !!!@ Isnt this redundant because of green?
+ !@ Isnt this redundant because of green?
INTEGER,INTENT(IN),OPTIONAL :: derive
!!!flag to toggle various derivatives of the solution (not to be used with
!!!green=ppm_poisson_grn_reprojec):
- !!!ppm_poisson_drv_none
- !!!ppm_poisson_drv_curl_fd (not implemented)
- !!!ppm_poisson_drv_grad_fd (not implemented)
- !!!ppm_poisson_drv_lapl_fd (not implemented)
- !!!ppm_poisson_drv_curl_sp (not implemented)
- !!!ppm_poisson_drv_grad_sp (not implemented)
- !!!ppm_poisson_drv_lapl_sp (not implemented)
+ !!! * ppm_poisson_drv_none
+ !!! * ppm_poisson_drv_curl_sp (not fully implemented)
+ !!! * ppm_poisson_drv_grad_sp (not implemented)
+ !!! * ppm_poisson_drv_lapl_sp (not implemented)
+ !!! * ppm_poisson_drv_div_sp (not implemented)
+ !!! * ppm_poisson_drv_curl_fd2
+ !!! * ppm_poisson_drv_grad_fd2 (not implemented)
+ !!! * ppm_poisson_drv_lapl_fd2 (not implemented)
+ !!! * ppm_poisson_drv_div_fd2 (not implemented)
+ !!! * ppm_poisson_drv_curl_fd4
+ !!! * ppm_poisson_drv_grad_fd4 (not implemented)
+ !!! * ppm_poisson_drv_lapl_fd4 (not implemented)
+ !!! * ppm_poisson_drv_div_fd4 (not implemented)
+ !!!
+ !!! curl=curl, grad=gradient,lapl=laplace operator,div=divergence
+ !!! sp=spectrally, fd2=2nd order finite differences, fd4=4th order FD
+ !!!
!!!The spectral derivatives can only be computed in this routine. Since
!!!the flag exists finite difference derivatives have also been included.
diff --git a/src/poisson/ppm_poisson_solve.f b/src/poisson/ppm_poisson_solve.f
index 585fc0f..5cbd1dc 100644
--- a/src/poisson/ppm_poisson_solve.f
+++ b/src/poisson/ppm_poisson_solve.f
@@ -1,35 +1,52 @@
!-------------------------------------------------------------------------
- ! ppm_poisson_solve.f90
+ ! Subroutine : ppm_poisson_solve.f90
+ !-------------------------------------------------------------------------
+ ! Copyright (c) 2010 CSE Lab (ETH Zurich), MOSAIC Group (ETH Zurich),
+ ! Center for Fluid Dynamics (DTU)
+ !
!-------------------------------------------------------------------------
!@mapping calls can be cleaned up
!-------------------------------------------------------------------------
-!!#define __ROUTINE ppm_poisson_solve
-!!#define __DIM 3
-!!#define __NCOM 3
-!!#define __ZEROSI (/0,0,0/)
-!#define __NOPE
SUBROUTINE __ROUTINE(topoid,meshid,ppmpoisson,fieldin,fieldout,gstw,info,&
& tmpcase)
+ !!! Routine to perform the Green's function solution of the Poisson
+ !!! equation. All settings are defined in ppm_poisson_initdef and stored
+ !!! in the ppmpoisson plan. The tmpcase argument allows the use of a
+ !!! different Green's function or operation than initialised. This is
+ !!! particularly useful for vorticity reprojection
+ !!! (ppm_poisson_grn_reprojec).
+ !!!
+ !!! The code will eventually get a needed cleanup overhaul.
USE ppm_module_map_field
USE ppm_module_map_field_global
USE ppm_module_map
-
+ USE ppm_module_typedef !@
+ USE ppm_module_data !@
+ USE ppm_module_finalize !@
IMPLICIT NONE
+ include 'mpif.h'
+
!-------------------------------------------------------------------------
! Arguments
!-------------------------------------------------------------------------
INTEGER, INTENT(IN) :: topoid
+ !!! Topology ID
INTEGER, INTENT(IN) :: meshid
+ !!! Mesh ID
TYPE(ppm_poisson_plan),INTENT(INOUT) :: ppmpoisson
- !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: fieldin
+ !!! The PPM Poisson plan
REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldin
- !REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER,INTENT(INOUT) :: fieldout
+ !!! Input data field
REAL(__PREC),DIMENSION(:,:,:,:,:),POINTER :: fieldout
+ !!! Output data field
INTEGER,DIMENSION(__DIM),INTENT(IN) :: gstw
+ !!! Ghost layer width
INTEGER, INTENT(OUT) :: info
+ !!! Return status, 0 upon succes
INTEGER,OPTIONAL,INTENT(IN) :: tmpcase
+ !!! Temporary operation (useful for ppm_poisson_grn_reprojec)
!-------------------------------------------------------------------------
! Local variables
@@ -45,6 +62,11 @@
REAL(__PREC) :: kx,ky,kz
REAL(__PREC) :: phix,phiy,phiz
REAL(__PREC) :: normfac
+ TYPE(ppm_t_equi_mesh), POINTER :: mesh => NULL()
+ TYPE(ppm_t_equi_mesh), POINTER :: target_mesh => NULL()
+ TYPE(ppm_t_topo), POINTER :: topo => NULL()
+ TYPE(ppm_t_topo), POINTER :: target_topo => NULL()
+
#ifndef __NOPE
INTEGER :: trank !@
trank =0
@@ -62,7 +84,6 @@ trank =0
ELSE
presentcase = ppmpoisson%case
ENDIF
-!@write(*,*) 'solving .... 1', ppm_rank
!-------------------------------------------------------------------------
!@ Perhaps check if ghostlayer suffices for a given fd stencil
@@ -70,88 +91,57 @@ trank =0
!-------------------------------------------------------------------------
- ! Perhaps allocate (and deallocate) arrays !@
+ ! Perhaps allocate (and deallocate) arrays
!-------------------------------------------------------------------------
-#ifdef __VARMESH
- !@tmp only works for one subdomain (not parallel)
- write(*,*) 'fieldin', ppm_rank
- DO isub=1,ppmpoisson%nsublistxy
- isubl=ppmpoisson%isublistxy(isub)
- DO k=1,ppmpoisson%ndataxy(3,isubl)/2-1
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataxy(1,isubl)/2-1
- DO j=1,ppmpoisson%ndataxy(2,isubl)/2-1
- write(*,'(A,E12.4,A,$)') ' (',fieldin(1,i,j,k,isub),')'
- END DO
- write(*,*)
- END DO
- END DO
- END DO
-#endif
- !-----------------------------------------------------------------------
- ! Set the real xy slabs 0 (part of the 0 padding) for free-space
- !@ free-space calculations and reprojection may cause problems !why?
- !-----------------------------------------------------------------------
- IF (presentcase .EQ. ppm_poisson_grn_pois_fre) THEN
- ppmpoisson%fldxyr = 0.0_MK
- ENDIF
- !-----------------------------------------------------------------------
- ! Map data globally to the slabs (XY)
- ! This is where the vorticity is extended and padded with 0 for free-space
- !-----------------------------------------------------------------------
- !Initialise
- CALL ppm_map_field_global(&
- & topoid, &
- & ppmpoisson%topoidxy, &
- & meshid, &
- & ppmpoisson%meshidxy,info)
- IF (info .NE. 0) THEN
- CALL ppm_write(ppm_rank,'ppm_poisson_solve','Failed to initialise field mapping.',info2)
- GOTO 9999
- ENDIF
- !Push the data
- CALL ppm_map_field_push(&
- & topoid, &
- & meshid,fieldin,__NCOM,info)
- IF (info .NE. 0) THEN
- CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to push vector field.',info2)
- GOTO 9999
- ENDIF
- !Send
- CALL ppm_map_field_send(info)
- IF (info .NE. 0) THEN
- CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to send field.',info2)
- GOTO 9999
- ENDIF
+ !-----------------------------------------------------------------------
+ ! Set the real xy slabs 0 (part of the 0 padding) for free-space
+ !@ free-space calculations and reprojection may cause problems !why?
+ !-----------------------------------------------------------------------
+ IF (presentcase .EQ. ppm_poisson_grn_pois_fre) THEN
+ ppmpoisson%fldxyr = 0.0_MK
+ ENDIF
+ !-----------------------------------------------------------------------
+ ! Map data globally to the slabs (XY)
+ ! This is where the vorticity is extended and padded with 0 for free-space
+ !-----------------------------------------------------------------------
+ !Initialise
+ CALL ppm_map_field_global(&
+ & topoid, &
+ & ppmpoisson%topoidxy, &
+ & meshid, &
+ & ppmpoisson%meshidxy,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank,'ppm_poisson_solve','Failed to initialise field mapping.',info2)
+ GOTO 9999
+ ENDIF
- !Retrieve
- CALL ppm_map_field_pop(&
- & ppmpoisson%topoidxy, &
- & ppmpoisson%meshidxy,ppmpoisson%fldxyr, &
- & __NCOM,__ZEROSI,info)
- IF (info .NE. 0) THEN
- CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
- GOTO 9999
- ENDIF
-!@write(*,*) 'solving .... 2', ppm_rank
-#ifdef __VARMESH
-!@tmp
- write(*,*) 'fldxyr', ppm_rank
- DO isub=1,ppmpoisson%nsublistxy
- isubl=ppmpoisson%isublistxy(isub)
- DO k=1,ppmpoisson%ndataxy(3,isubl)-1
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataxy(1,isubl)-1
- DO j=1,ppmpoisson%ndataxy(2,isubl)-1
- write(*,'(A,E12.4,A,$)') ' (',ppmpoisson%fldxyr(1,i,j,k,isub),')'
- END DO
- write(*,*)
- END DO
- END DO
- END DO
-#endif
+ !Push the data
+ CALL ppm_map_field_push(&
+ & topoid, &
+ & meshid,fieldin,__NCOM,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to push vector field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Send
+ CALL ppm_map_field_send(info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to send field.',info2)
+ GOTO 9999
+ ENDIF
+
+ !Retrieve
+ CALL ppm_map_field_pop(&
+ & ppmpoisson%topoidxy, &
+ & ppmpoisson%meshidxy,ppmpoisson%fldxyr, &
+ & __NCOM,__ZEROSI,info)
+ IF (info .NE. 0) THEN
+ CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
+ GOTO 9999
+ ENDIF
!-----------------------------------------------------------------------
! Do slab FFT (XY) - use the non-reduced topology !@what does this mean
@@ -160,7 +150,6 @@ trank =0
& ppmpoisson%meshidxy, ppmpoisson%planfxy, &
& ppmpoisson%fldxyr, ppmpoisson%fldxyc, &
& info)
-!@write(*,*) 'solving .... 3', ppm_rank
#ifdef __NOPE
if (ppm_rank .EQ. trank) THEN
@@ -250,7 +239,6 @@ trank =0
END DO
ENDIF
#endif
-!@write(*,*) 'solving .... 4', ppm_rank
!-----------------------------------------------------------------------
! Map to the pencils (Z)
@@ -293,7 +281,6 @@ trank =0
CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
GOTO 9999
ENDIF
-!@write(*,*) 'solving .... 5', ppm_rank
!-----------------------------------------------------------------------
! Do pencil FFT (Z)
@@ -303,7 +290,6 @@ trank =0
& ppmpoisson%fldzc1, ppmpoisson%fldzc2, &
& info)
-!@write(*,*) 'solving .... 6', ppm_rank
#ifdef __NOPE
if (ppm_rank .EQ. trank) THEN
write(*,*)
@@ -365,7 +351,6 @@ trank =0
END DO
ENDIF
#endif
-!@write(*,*) 'solving .... 7', ppm_rank
!-----------------------------------------------------------------------
! Apply the periodic Green's function
@@ -401,14 +386,10 @@ trank =0
& ppmpoisson%fldzc2(2,i,j,k,isub)
ppmpoisson%fldzc2(3,i,j,k,isub) = ppmpoisson%fldgrnc( i,j,k,isub)*&
& ppmpoisson%fldzc2(3,i,j,k,isub)
- !!ppmpoisson%fldzc2(1,i,j,k,isub) = ppmpoisson%fldgrnc( i,j,k,isub)
- !!ppmpoisson%fldzc2(2,i,j,k,isub) = ppmpoisson%fldgrnc( i,j,k,isub)
- !!ppmpoisson%fldzc2(3,i,j,k,isub) = ppmpoisson%fldgrnc( i,j,k,isub)
ENDDO
ENDDO
ENDDO
ENDDO
-!@write(*,*) 'solving .... 8', ppm_rank
!-----------------------------------------------------------------------
! Spectral derivatives
! normkx, etc contains 2pi/Lx
@@ -452,14 +433,12 @@ trank =0
!-----------------------------------------------------------------------
! Vorticity re-projection
- !@ has the domain length really been included in these wave numbers
!-----------------------------------------------------------------------
ELSE IF (presentcase .EQ. ppm_poisson_grn_reprojec) THEN
!remembering to normalize the FFT
normfac = 1.0_MK/ REAL((ppmpoisson%nmz(1)-1)* & !vertex
& (ppmpoisson%nmz(2)-1)* &
& (ppmpoisson%nmz(3)-1),MK)
- write(*,*) 'reprojection ppm style' !@
DO isub=1,ppmpoisson%nsublistz
isubl=ppmpoisson%isublistz(isub)
DO k=1,ppmpoisson%ndataz(3,isubl)
@@ -496,7 +475,6 @@ trank =0
ENDDO
ENDIF
-!@write(*,*) 'solving .... 9', ppm_rank
!-----------------------------------------------------------------------
! IFFT pencil (Z)
!-----------------------------------------------------------------------
@@ -505,7 +483,6 @@ trank =0
& ppmpoisson%fldzc2, ppmpoisson%fldzc1, &
& info)
-!@write(*,*) 'solving .... 11', ppm_rank
#ifdef __NOPE
if (ppm_rank .EQ. trank) THEN
write(*,*)
@@ -592,7 +569,6 @@ trank =0
GOTO 9999
ENDIF
-!@write(*,*) 'solving .... 12', ppm_rank
!-----------------------------------------------------------------------
! IFFT (XY) use the non-reduced topology
!-----------------------------------------------------------------------
@@ -601,7 +577,6 @@ trank =0
& ppmpoisson%fldxyc, ppmpoisson%fldxyr, &
& info)
-!@write(*,*) 'solving .... 13', ppm_rank
#ifdef __NOPE
if (ppm_rank .EQ. trank) THEN
write(*,*)
@@ -648,22 +623,6 @@ trank =0
ENDIF
#endif
-#ifdef __VARMESH
-!@tmp
- write(*,*) 'fldxyr2', ppm_rank
- DO isub=1,ppmpoisson%nsublistxy
- isubl=ppmpoisson%isublistxy(isub)
- DO k=1,ppmpoisson%ndataxy(3,isubl)
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataxy(1,isubl)
- DO j=1,ppmpoisson%ndataxy(2,isubl)
- write(*,'(A,E12.4,A,$)') ' (',ppmpoisson%fldxyr(1,i,j,k,isub),')'
- END DO
- write(*,*)
- END DO
- END DO
- END DO
-#endif
!-----------------------------------------------------------------------
! Map back to standard topology (XYZ)
!-----------------------------------------------------------------------
@@ -677,7 +636,7 @@ trank =0
CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to initialise field mapping.',info2)
GOTO 9999
ENDIF
-!@write(*,*) 'solving .... 13a', ppm_rank
+
!Push the data
CALL ppm_map_field_push(&
@@ -687,7 +646,20 @@ trank =0
CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to push vector field.',info2)
GOTO 9999
ENDIF
-!@write(*,*) 'solving .... 13b', ppm_rank, ppmpoisson%topoidxy, topoid, ppmpoisson%meshidxy, meshid
+ topo => ppm_topo(topoid)%t!@
+ mesh => topo%mesh(meshid) !@
+ target_topo => ppm_topo(ppmpoisson%topoidxy)%t !@
+ target_mesh => target_topo%mesh(ppmpoisson%meshidxy) !@
+ DO isub=1,topo%nsublist
+ isubl=topo%isublist(isub)
+ !@write(*,*) 'johannestest rank', ppm_rank,'istart', mesh%istart(:,isubl), &
+ !@'nnodes',mesh%nnodes(:,isubl), mesh%nm
+ ENDDO
+ DO isub=1,ppmpoisson%nsublistxy
+ isubl=ppmpoisson%isublistxy(isub)
+ !@write(*,*) 'johannestestxy rank', ppm_rank,'istart', target_mesh%istart(:,isubl), &
+ !@'nnodes', target_mesh%nnodes(:,isubl), target_mesh%nm
+ ENDDO
!Send
CALL ppm_map_field_send(info)
@@ -696,7 +668,6 @@ trank =0
GOTO 9999
ENDIF
-!@write(*,*) 'solving .... 14', ppm_rank
!-------------------------------------------------------------------------
! FINAL RETRIEVE - Here we do different things depending on the task
! i.e. the receiver varies
@@ -705,38 +676,32 @@ trank =0
ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_fd4) .AND. &
& (presentcase .EQ. ppm_poisson_grn_pois_per .OR. & !@these may be unnecessary - perhaps just the derive value. Or maybe not: in case of vorticity reprojection we could get lost
presentcase .EQ. ppm_poisson_grn_pois_fre)) THEN
-!@write(*,*) 'solving .... 14a', ppm_rank
CALL ppm_map_field_pop(&
& topoid, &
& meshid,ppmpoisson%drv_vr, &
- & __NCOM,gstw,info) !!! gst
+ & __NCOM,gstw,info)
IF (info .NE. 0) THEN
CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
GOTO 9999
ENDIF
-!@write(*,*) 'solving .... 14b', ppm_rank
!-------------------------------------------------------------------------
! Ghost the temporary array for derivatives (drv_vr)
!-------------------------------------------------------------------------
-!@write(*,*) 'solving .... 14c', ppm_rank
CALL ppm_map_field_ghost_get(topoid,meshid,gstw,info)
IF (info .NE. 0) THEN
CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to initialise ghosts.',info2)
GOTO 9999
ENDIF
-!@write(*,*) 'solving .... 14d', ppm_rank
CALL ppm_map_field_push(topoid,meshid,ppmpoisson%drv_vr,__NCOM,info)
IF (info .NE. 0) THEN
CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to push ghosts.',info2)
GOTO 9999
ENDIF
-!@write(*,*) 'solving .... 14e', ppm_rank
CALL ppm_map_field_send(info)
IF (info .NE. 0) THEN
CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to send ghosts.',info2)
GOTO 9999
ENDIF
-!@write(*,*) 'solving .... 14f', ppm_rank
CALL ppm_map_field_pop(topoid,meshid,ppmpoisson%drv_vr,__NCOM,gstw,info)
IF (info .NE. 0) THEN
CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop ghosts.',info2)
@@ -744,38 +709,19 @@ trank =0
ENDIF
ELSE
-!@write(*,*) 'solving .... 14g', ppm_rank
CALL ppm_map_field_pop(&
& topoid, &
& meshid,fieldout, &
- & __NCOM,gstw,info) !!! gst
-!@write(*,*) 'solving .... 14h', ppm_rank
+ & __NCOM,gstw,info)
ENDIF
IF (info .NE. 0) THEN
CALL ppm_write(ppm_rank, 'ppm_poisson_solve','Failed to pop vector field.',info2)
GOTO 9999
ENDIF
-!@write(*,*) 'solving .... 15', ppm_rank
-#ifdef __VARMESH
-!@tmp only works for one subdomain (not parallel)
- write(*,*) 'fieldout', ppm_rank
- DO isub=1,ppmpoisson%nsublistxy
- isubl=ppmpoisson%isublistxy(isub)
- DO k=1,ppmpoisson%ndataxy(3,isubl)/2-1
- write(*,*) 'z',k
- DO i=1,ppmpoisson%ndataxy(1,isubl)/2-1
- DO j=1,ppmpoisson%ndataxy(2,isubl)/2-1
- write(*,'(A,E12.4,A,$)') ' (',fieldout(1,i,j,k,isub),')'
- END DO
- write(*,*)
- END DO
- END DO
- END DO
-#endif
!-------------------------------------------------------------------------
- !@To come: treat ghost layer to make FD stencils work
+ ! Treat ghost layer to make FD stencils work
!-------------------------------------------------------------------------
IF (ppmpoisson%derivatives .EQ. ppm_poisson_drv_curl_fd2 .AND.&
& (presentcase .EQ. ppm_poisson_grn_pois_fre)) THEN
@@ -788,7 +734,6 @@ trank =0
& 2,4,gstw,info)
ENDIF
-!@write(*,*) 'solving .... 16', ppm_rank
!-------------------------------------------------------------------------
! Optionally do derivatives
! Perhaps make ppm_poisson_fd take _none as argument. Then maybe no
@@ -807,7 +752,6 @@ trank =0
& ppm_poisson_drv_curl_fd4,info)
ENDIF
-!@write(*,*) 'solving .... 17', ppm_rank
!-------------------------------------------------------------------------
! Finally ghost the velocity/stream function field before returning it
! Also extrapolate if freespace
@@ -821,7 +765,6 @@ trank =0
& 2,4,gstw,info)
ENDIF
-!@write(*,*) 'solving .... 18', ppm_rank
!-------------------------------------------------------------------------
! Perhaps allocate (and deallocate) arrays !@
!-------------------------------------------------------------------------
diff --git a/src/ppm_define.h b/src/ppm_define.h
index 91803f6..9526c22 100644
--- a/src/ppm_define.h
+++ b/src/ppm_define.h
@@ -1,2 +1,3 @@
#define __MPI
#define __FFTW
+#define __Linux
diff --git a/src/ppm_module_poisson.f b/src/ppm_module_poisson.f
index 44e0f35..fd0408d 100644
--- a/src/ppm_module_poisson.f
+++ b/src/ppm_module_poisson.f
@@ -18,9 +18,13 @@
! derivatives (optional): none,curl,gradient,laplace, - spectral/FD
! possibly the option to de- and reallocate arrays of slabs/pencils
!
+ ! The fields returned from the subroutine have been ghosted/extrapolated
+ ! as necessary.
!
! Add finalize routine
! This routine needs to be modified to allow cell centred data
+ ! These routines should all be renamed - potentially not just Poisson
+ ! equations can be solved.
!-------------------------------------------------------------------------
#define __SINGLE 0
#define __DOUBLE 1
--
GitLab