Skip to content
Snippets Groups Projects
Commit 83ac556d authored by jtra's avatar jtra
Browse files

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
parent 84c9c96c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/poisson/ppm_poisson_extrapolateghost.f
+ 20
− 12
View file @ 83ac556d
!-------------------------------------------------------------------------
! 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
!-------------------------------------------------------------------------
......
src/poisson/ppm_poisson_fd.f
+ 35
− 8
View file @ 83ac556d
!-------------------------------------------------------------------------
! 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
!-------------------------------------------------------------------------
......
src/poisson/ppm_poisson_init_predef.f
+ 47
− 18
View file @ 83ac556d
!-------------------------------------------------------------------------
! 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.
......
src/poisson/ppm_poisson_solve.f
+ 96
− 153
View file @ 83ac556d
!-------------------------------------------------------------------------
! 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 !@
!-------------------------------------------------------------------------
......
src/ppm_define.h
+ 1
− 0
View file @ 83ac556d
#define __MPI
#define __FFTW
#define __Linux
src/ppm_module_poisson.f
+ 4
− 0
View file @ 83ac556d
......@@ -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
......
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