diff --git a/HySoP/CMakeLists.txt b/HySoP/CMakeLists.txt
index 71c2db0426cbbced648c67da0bb84d712f172a9c..3d665d59ddd1d4609b3252e2e4217debe40ac33e 100644
--- a/HySoP/CMakeLists.txt
+++ b/HySoP/CMakeLists.txt
@@ -28,7 +28,7 @@ option(VERBOSE_MODE "enable verbose mode for cmake exec. Default = on" ON)
option(USE_PPM "link parmes with ppm when this mode is enable. Default = on." ON)
option(USE_MPI "compile and link parmes with mpi when this mode is enable. Default = on." ON)
option(USE_PPM "link Parmes with PPM library (core component). Default = on." ON)
-option(USE_PPMNumerics "link Parmes with PPM-numerics. Default = off" OFF)
+option(USE_PPMNumerics "link Parmes with PPM-numerics. Default = off" ON)
option(DOUBLEPREC "set precision for real numbers to double precision when this mode is enable. Default = on." ON)
option(WITH_TESTS "Enable testing. Default = on" ON)
option(BUILD_SHARED_LIBS "Enable dynamic library build, default = ON" ON)
diff --git a/HySoP/src/interfaces/ppm/wrap_ppm_topologies.f95 b/HySoP/src/interfaces/ppm/wrap_ppm_topologies.f95
index 3f9ee30247d00dae9584cdb5510ca1a264bc328e..a13e7dc1f34cd0d92762a7a64108f51a825f9472 100644
--- a/HySoP/src/interfaces/ppm/wrap_ppm_topologies.f95
+++ b/HySoP/src/interfaces/ppm/wrap_ppm_topologies.f95
@@ -19,13 +19,13 @@ contains
type(c_Ptr), intent(in), VALUE :: minPhys
type(c_Ptr), intent(in), VALUE :: maxPhys
type(c_Ptr), intent(in), VALUE :: bc
- real(kind=real_kind), intent(in) :: ghostsize
+ real(kind=real_kind), intent(in) :: ghostsize ! ghostsize is a real in ppm subroutine ... strange ...
! type(c_ptr), intent(inout) :: costPerProc
! Local vars
integer :: info
- integer :: nbProcs = 3 ! TODO : input arg
+ ! integer :: nbProcs = 3 ! TODO : input arg
integer :: assig = ppm_param_assign_internal
real(kind=real_kind), pointer, dimension(:) :: cost => NULL()
real(kind=real_kind), pointer, dimension(:) ::min_phys => NULL(), max_phys => NULL()
diff --git a/HySoP/src/main/Domain.f90 b/HySoP/src/main/Domain.f90
index e84fb57ebc8aa6bdeabb14a9e00f6659a67d1954..039f9e3f498ec801e49af2a6f13498ac615564ba 100755
--- a/HySoP/src/main/Domain.f90
+++ b/HySoP/src/main/Domain.f90
@@ -11,8 +11,10 @@ module Domain
integer, dimension(:), pointer :: domain_ghostsize=>NULL()
integer, dimension(:), pointer :: grid_resolution =>NULL()
+ real(mk), dimension(:), pointer :: grid_step =>NULL()
integer, private :: istat
+
contains
subroutine init_geometry()
@@ -21,7 +23,7 @@ contains
allocate(domain_minCoords(dime), domain_maxCoords(dime), stat = istat)
if(istat.ne.0) stop 'Geometry, allocation failed'
domain_minCoords = 0.0
- domain_maxCoords = 3.1415927
+ domain_maxCoords = 1.1
! Boundary conditions and ghosts
allocate(domain_bc(2*dime), domain_ghostsize(dime), stat = istat)
if(istat.ne.0) stop 'BC, allocation failed'
@@ -32,9 +34,10 @@ contains
subroutine init_grid()
- allocate(grid_resolution(dime),stat=istat)
+ allocate(grid_resolution(dime),grid_step(dime),stat=istat)
if(istat.ne.0) stop 'grid_resolution, allocation failed'
- grid_resolution = 10
+ grid_resolution = 65
+ grid_step = (domain_maxCoords - domain_minCoords)/(real(grid_resolution,mk)-1.)
end subroutine init_grid
diff --git a/HySoP/src/main/Fields.f90 b/HySoP/src/main/Fields.f90
index d0e7a4a705cc87aa17fa2004cdcc81107c6cb5d5..679adad9b960263f40dca70f53175189e26b34d6 100755
--- a/HySoP/src/main/Fields.f90
+++ b/HySoP/src/main/Fields.f90
@@ -2,22 +2,23 @@ module Fields
use client_data, only: mk,dime
use ppm_module_typedef, only : ppm_t_equi_mesh, ppm_t_topo
-
+
implicit none
!> Velocity
real(mk), dimension(:,:,:,:,:), pointer :: velocity => NULL()
!> Vorticity
real(mk), dimension(:,:,:,:,:), pointer :: vorticity =>NULL()
+ !> Stream function
+ real(mk), dimension(:,:,:,:,:), pointer :: stream_function =>NULL()
contains
- subroutine init_fields(ghostsize, mesh, topo)
+ subroutine init_fields(ghostsize, topo)
integer, dimension(:), pointer :: ghostsize
- type(ppm_t_equi_mesh), pointer :: mesh
type(ppm_t_topo), pointer :: topo
-
+
!> Number of points in each dir (1st index) for each sub (2nd index)
integer, dimension(:,:), pointer :: ndata =>NULL()
!> List of subdomains for the current proc - Warning, it seems that from ppm, isublist is of size nbproc (total number of proc) but
@@ -31,11 +32,11 @@ contains
! Lower and upper bounds for fields
integer, dimension(dime) :: ldl, ldu
integer, dimension(dime) :: vndata
-
- ndata => mesh%nnodes
+
+ ndata => topo%mesh(1)%nnodes
isublist => topo%isublist
nsublist = topo%nsublist
-
+
vndata(1) = maxval(ndata(1,isublist(1:nsublist)))
vndata(2) = maxval(ndata(2,isublist(1:nsublist)))
vndata(3) = maxval(ndata(3,isublist(1:nsublist)))
@@ -48,40 +49,12 @@ contains
! Vorticity
allocate(vorticity(dime,ldl(1):ldu(1),ldl(2):ldu(2),ldl(3):ldu(3),nsublist), stat = istat)
if(istat.ne.0) stop 'Field allocation error for vorticity'
+ ! Stream function
+ allocate(stream_function(dime,ldl(1):ldu(1),ldl(2):ldu(2),ldl(3):ldu(3),nsublist), stat = istat)
+ if(istat.ne.0) stop 'Field allocation error for stream_function'
- print *, shape(velocity)
nullify(isublist, ndata)
end subroutine init_fields
-
-!!$ subroutine allocfield(field, field_length, info)
-!!$
-!!$ real(mk), dimension(:,:,:,:,:), pointer :: field
-!!$ integer, intent(out) :: info
-!!$ integer, intent(in ) :: field_length
-!!$
-!!$ integer, dimension(3) :: vndata, ldl, ldu
-!!$
-!!$ info = -1
-!!$
-!!$ ! Check if field is already allocated and nullify if required.
-!!$ if(associated(field)) then
-!!$ deallocate(field,stat=info)
-!!$ nullify(field)
-!!$ if(info.ne.0) stop 'tools:allocfield allocation failed'
-!!$ end if
-!!$
-!!$
-!!$ vndata(1) = maxval(ndata(1,isublist(1:nsublist)))
-!!$ vndata(2) = maxval(ndata(2,isublist(1:nsublist)))
-!!$ vndata(3) = maxval(ndata(3,isublist(1:nsublist)))
-!!$
-!!$ ldl = 1 - ghostsize
-!!$ ldu = vndata + ghostsize
-!!$
-!!$ allocate(field(field_length,ldl(1):ldu(1),ldl(2):ldu(2),ldl(3):ldu(3),nsublist),&
-!!$ & stat = info)
-!!$
-!! end subroutine allocfield
-
+
end module Fields
diff --git a/HySoP/src/main/FieldsComputation.f90 b/HySoP/src/main/FieldsComputation.f90
new file mode 100755
index 0000000000000000000000000000000000000000..c4f415bc5c71397960c7a9b2c2e1f1dfb8dd8b2c
--- /dev/null
+++ b/HySoP/src/main/FieldsComputation.f90
@@ -0,0 +1,42 @@
+module FieldsComputation
+
+ use client_data, only: mk
+
+ !use Solver, only: solve_poisson
+
+ implicit none
+
+ private
+
+ public :: fit_velocity
+
+contains
+
+!!$ subroutine compute_velocity(field_out, field_in,topoid,meshid, ghostsize)
+!!$ real(mk), dimension(:,:,:,:,:), pointer :: field_out
+!!$ !> Vorticity
+!!$ real(mk), dimension(:,:,:,:,:), pointer :: field_in
+!!$ integer, intent(in) :: topoid, meshid
+!!$
+!!$ integer, dimension(:), pointer :: ghostsize
+!!$
+!!$ ! Solves poisson and computes velocity
+!!$ call solve_poisson(field_in, field_out, topoid, meshid,ghostsize)
+!!$
+!!$ ! Multigrid in PPM is bugged ...
+!!$ !call solve_poisson(topoid, stream_function, field_in)
+!!$
+!!$ end subroutine compute_velocity
+
+ subroutine fit_velocity(vel, resolution, bc)
+ ! Velocity, intent(inout)
+ real(mk), dimension(:,:,:,:,:), pointer :: vel
+ ! grid resolution, intent(in)
+ integer, dimension(:), pointer :: resolution
+ ! boundary conditions type, intent(in)
+ integer, dimension(:), pointer :: bc
+
+ ! TODO
+ end subroutine fit_velocity
+
+end module FieldsComputation
diff --git a/HySoP/src/main/NavierStokes3D.f90 b/HySoP/src/main/NavierStokes3D.f90
new file mode 100755
index 0000000000000000000000000000000000000000..54c56eb77db902960159a15bf0535d1198937338
--- /dev/null
+++ b/HySoP/src/main/NavierStokes3D.f90
@@ -0,0 +1,158 @@
+!> Temp modules for ppm_client
+module ppmNavierStokes3D
+
+ use ppm_module_init
+ use ppm_module_data, only : ppm_kind_double
+ use ppm_module_finalize
+
+ ! use client_io
+ use client_data, only: mk, dime
+ ! Physical domain and grid
+ use Domain
+ ! Fields on the grid
+ use Fields, only: init_fields, velocity, vorticity
+ ! Topology
+ use client_topology, only: init_topo, mesh, topo
+ ! Multigrid solver
+ use Solver, only : init_poisson_solver, solve_poisson
+
+ use FieldsComputation, only : fit_velocity
+
+ use mpi
+ !use WrapFort
+
+ implicit none
+
+ include "ppm_numerics.h"
+
+ integer, private :: info
+
+contains
+
+ !> All required initialisations :
+ !! MPI, ppm
+ !! Domain, grid
+ !! Fields
+ !! Particles
+ subroutine init_client()
+
+ integer :: prec,tol
+ ! MPI comm
+ integer :: comm
+ ! debug mode
+ integer :: debug
+ ! error status
+ integer :: info
+
+ !======================
+ ! Init ppm
+ !======================
+ prec = ppm_kind_double ! Defined in ppm_param.h
+ comm = MPI_COMM_WORLD !
+ debug = 2
+ tol = -10
+ info = -1
+ call ppm_init(dime,prec,tol,comm,debug,info)
+
+ !======================
+ ! Geometry and grid
+ !======================
+ call init_geometry()
+ call init_grid()
+
+ !======================
+ ! Creates the topology
+ !======================
+ call init_topo(domain_minCoords, domain_maxCoords, domain_bc, domain_ghostsize, grid_resolution)
+
+ !======================
+ ! Fields allocation
+ !======================
+ call init_fields(domain_ghostsize, topo)
+
+ !======================
+ ! Init Physics
+ !======================
+ velocity = 0.0_mk
+ vorticity = 1.0_mk
+
+ vorticity(3,10,10,10,:) = 50.0_mk
+
+ !======================
+ ! Init Particles
+ !======================
+
+ print *, "end of parmes:ppm:init_client"
+
+ end subroutine init_client
+
+ subroutine main_client() bind(c,name='NavierStokes3D')
+
+ real(mk) :: initial_time, final_time, current_time, time_step
+ integer :: i, rank
+ initial_time = 0.0
+ final_time = 1.0
+ time_step = 0.005
+ current_time = initial_time
+
+ call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info)
+
+ ! Initialisation stuff
+ call init_client()
+
+ !======================
+ ! Init solver
+ !======================
+ ! Init multigrid works in ppm (or seems to ...) but solve fails.
+ !call init_multigrid(topo%ID, topo%mesh(1)%ID, domain_ghostsize, domain_bc)
+
+ call init_poisson_solver(vorticity,velocity,topo%ID, topo%mesh(1)%ID)
+
+!!$ print *, shape(topo%mesh)
+!!$ do i=1, size(topo%mesh)
+!!$ print *, 'mesh :', topo%mesh(i)%ID
+!!$ end do
+!!$
+!!$ do i=1, 3!size(ppm_topo)
+!!$ print *, rank, ' topo: ', ppm_topo(i)%t%ID
+!!$
+!!$ print *, rank, ' minmax: ', ppm_topo(i)%t%min_subd, ' max ...', ppm_topo(i)%t%max_subd
+!!$ print *, rank ,' bc ', ppm_topo(i)%t%subs_bc
+!!$
+!!$ end do
+!!$
+!!$
+ ! Time loop
+
+ do while(current_time < final_time )
+
+ ! Compute velocity from vorticity
+ ! Two steps:
+ ! - solve Poisson for stream_function and update velocity from stream_function.
+ ! - update velocity to fit with required boundary values
+ ! call solve_poisson(vorticity, velocity, topo%ID, topo%mesh(1)%ID,domain_ghostsize)
+ ! call fit_velocity(velocity, grid_resolution,domain_bc)
+
+ ! Penalize
+! call penalize_velocity(velocity)
+
+ ! Compute the new "penalized" vorticity
+ !call compute_vorticity(vorticity, velocity)
+
+ ! Initialize the particles
+ !call create_particles(vorticity,velocity,wp,xp,vp)
+
+ ! Integrate
+ !call push_particles(xp,vp,wp)
+ ! Remesh
+ !call remesh_particles(wp,vorticity)
+
+ current_time = current_time + time_step
+ current_time = final_time
+ enddo
+
+ call ppm_finalize(info)
+ print *, 'end ppm simulation'
+ end subroutine main_client
+
+end module ppmNavierStokes3D
diff --git a/HySoP/src/main/Penalize.f90 b/HySoP/src/main/Penalize.f90
new file mode 100755
index 0000000000000000000000000000000000000000..2b9b1856cb54142fe1fe5d1a1605a6f5d7ccc403
--- /dev/null
+++ b/HySoP/src/main/Penalize.f90
@@ -0,0 +1,124 @@
+module penalisation
+
+ use client_data, only: mk,dime
+ use ppm_module_typedef, only : ppm_t_equi_mesh, ppm_t_topo
+
+ implicit none
+
+ private
+
+ public :: penalisation_init, penalise_velocity
+
+ ! Penalisation kernel function
+ real(mk), dimension(:,:,:,:), pointer :: chi => NULL()
+ !
+ real(mk) :: lambda
+
+contains
+
+ subroutine penalisation_init(topo, ghostsize, meshStep)
+ integer, dimension(:), pointer :: ghostsize
+ type(ppm_t_topo), pointer :: topo
+ integer, dimension(:), pointer :: meshStep
+
+ !> Number of points in each dir (1st index) for each sub (2nd index)
+ integer, dimension(:,:), pointer :: ndata =>NULL()
+ !> List of subdomains for the current proc - Warning, it seems that from ppm, isublist is of size nbproc (total number of proc) but
+ ! for each proc, only isublist(1:nsublist) is relevent
+ integer, dimension(:),pointer :: isublist => NULL()
+ !> Number of subs for the current proc
+ integer :: nsublist
+
+ integer::istat,i,j,k,isub,isubl
+
+ ! Lower and upper bounds for fields
+ integer, dimension(dime) :: ldl, ldu
+ integer, dimension(dime) :: vndata
+
+ ! z position
+ real(mk) :: zPos
+ ! lower bound of the physical domain
+ real(mk), dimension(:), pointer :: minPhys => NULL()
+ ! mesh indexes
+ integer, dimension(:,:), pointer :: istart => NULL()
+ ! virtual boundary layer position
+ real(mk) :: limitPos
+ limitPos = 1.0
+
+ ! Set lambda value
+ lambda = 5.5D5
+
+ ndata => topo%mesh(1)%nnodes
+ isublist => topo%isublist
+ nsublist = topo%nsublist
+
+ vndata(1) = maxval(ndata(1,isublist(1:nsublist)))
+ vndata(2) = maxval(ndata(2,isublist(1:nsublist)))
+ vndata(3) = maxval(ndata(3,isublist(1:nsublist)))
+
+ ldl = 1 - ghostsize
+ ldu = vndata + ghostsize
+
+ allocate(chi(ldl(1):ldu(1),ldl(2):ldu(2),ldl(3):ldu(3),nsublist),stat = istat)
+ if(istat.ne.0) stop 'Chi function allocation error.'
+
+ chi = 0.0
+ istart => topo%mesh(1)%istart
+ minPhys => topo%min_physd
+ ! Loop over subprocess ...
+ do isub=1,nsublist
+ ! loop over mesh nodes
+ isubl=isublist(isub)
+ do k=1,ndata(3,isubl)
+ do j=1,ndata(2,isubl)
+ do i=1,ndata(1,isubl)
+ zpos = minPhys(3) + meshStep(3)*istart(k,isubl)
+ if( (zpos>limitPos).or.(zpos<-limitPos)) then
+ chi(i,j,k,isub) = lambda
+ end if
+ end do
+ end do
+ end do
+ end do
+
+ nullify(isublist, ndata)
+
+
+ end subroutine penalisation_init
+ !>
+ subroutine penalise_velocity(vel, topo, resolution, dt)
+ ! Velocity, intent(inout)
+ real(mk), dimension(:,:,:,:,:), pointer :: vel
+ ! The topology
+ type(ppm_t_topo), pointer :: topo
+ ! grid resolution, intent(in)
+ integer, dimension(:), pointer :: resolution
+ ! Time step
+ real(mk), intent(in) :: dt
+ ! Chi penalisation function
+
+
+ integer :: i,j,k, isub, nsublist,isubl
+ real(mk) :: coef
+
+ nsublist = topo%nsublist
+
+
+ do isub = 1, nsublist
+ isubl=topo%isublist(isub)
+ do k=1,resolution(3)
+ do j=1,resolution(2)
+ do i=1,resolution(1)
+
+ coef =1./(1.+ dt*chi(i,j,k,isubl))
+ vel(1,i,j,k,isubl)=vel(1,i,j,k,isubl)*coef
+ vel(2,i,j,k,isubl)=vel(2,i,j,k,isubl)*coef
+ vel(3,i,j,k,isubl)=vel(3,i,j,k,isubl)*coef
+
+ end do
+ end do
+ end do
+ end do
+end subroutine penalise_velocity
+
+end module Penalisation
diff --git a/HySoP/src/main/Plouhmans.f90 b/HySoP/src/main/Plouhmans.f90
index 7127ac1ddcefe474c1fd99a0d3cda1a9aee9a91d..6469e385f148628459c411f5dbf3356e4f86e54f 100644
--- a/HySoP/src/main/Plouhmans.f90
+++ b/HySoP/src/main/Plouhmans.f90
@@ -4,21 +4,20 @@ module ppmExample
use ppm_module_init
use ppm_module_data, only : ppm_kind_double
use ppm_module_finalize
- use ppm_module_mg_init
- ! For ppm_param
-! use client_io
+
+ ! use client_io
use client_data, only: mk, dime
-
! Physical domain and grid
use Domain
- use Fields, only: init_fields
+ ! Fields on the grid
+ use Fields, only: init_fields, velocity, vorticity, stream_function
! Topology
use client_topology, only: init_topo, mesh, topo
+ ! Multigrid solver
+ !use Solver, only : init_multigrid, solve_poisson
- use Solver, only : init_multigrid
use mpi
use WrapFort
- use Tools
implicit none
@@ -58,6 +57,7 @@ contains
!======================
call init_geometry()
call init_grid()
+
!======================
! Creates the topology
!======================
@@ -66,23 +66,36 @@ contains
!======================
! Fields allocation
!======================
- call init_fields(domain_ghostsize, mesh, topo)
+ call init_fields(domain_ghostsize, topo)
!======================
! Init solver
!======================
- call init_multigrid(topo%ID, mesh%ID, domain_ghostsize, domain_bc)
+! call init_multigrid(topo%ID, mesh%ID, domain_ghostsize, domain_bc)
+ !======================
+ ! Init Physics
+ !======================
+ velocity = 0.0_mk
+ vorticity = 0.0_mk
+
+ !======================
+ ! Init Particles
+ !======================
+
print *, "end of parmes:ppm:init_client"
end subroutine init_client
subroutine main_client() bind(c,name='plouhmans')
+ ! Multigrid parameters ...
print *, 'run ppm simulation ...'
! init ppm ...
call init_client()
+! call solve_poisson(topo%ID, stream_function, vorticity)
+
call ppm_finalize(info)
print *, 'end ppm simulation'
end subroutine main_client
diff --git a/HySoP/src/main/Solver.f90 b/HySoP/src/main/Solver.f90
index 9a23eca298ad27ef49e1ce021a575092c491bbd8..70a7582d3df807c110bedd94dfa804bfebabbb22 100755
--- a/HySoP/src/main/Solver.f90
+++ b/HySoP/src/main/Solver.f90
@@ -4,9 +4,32 @@ module Solver
use client_data, only: mk, dime
use ppm_module_mg_init
use ppm_module_numerics_data, only: ppm_param_eq_poisson, ppm_param_smooth_rbsor
+ !use ppm_module_fdsolver_init
+ !use ppm_module_fdsolver_solve
+ use ppm_module_mg_solv
+ use ppm_module_poisson, only: ppm_poisson_init, ppm_poisson_plan, ppm_poisson_grn_pois_per, ppm_poisson_solve, &
+ ppm_poisson_drv_curl_fd2, ppm_poisson_fd
+ use client_topology
+ use ppm_module_typedef, only : ppm_param_mesh_coarsen
+ use ppm_module_mesh_derive, only: ppm_mesh_derive
implicit none
+ private
+
+ public :: init_poisson_solver, solve_poisson
+
+ type(ppm_poisson_plan), pointer :: fftwplanForppm => NULL()
+
+ ! Interface to poisson solver initialisation. Only fft at the time
+ interface init_poisson_solver
+ module procedure init_fftw
+ end interface
+
+ interface solve_poisson
+ module procedure solve_poisson_fftw
+ end interface
+
contains
subroutine init_multigrid(topoid, meshid, ghostsize, bc)
@@ -22,23 +45,12 @@ contains
logical :: with_w_cycles, dump_info
real(mk) :: omega
integer :: limlev ! Number of mg levels
- integer :: info, i,j,k,l
- CHARACTER(LEN=120) :: mesg
- i =1
- j=2
- k=3
- l=4
-
- WRITE(mesg,'(A,I5,A,I1,A,I6,A,I3)') 'Mesh on sub ', &
- i,' in dimension ',j,' has ',k, &
- ' mesh cells and is not divisible by ',l
-
- print *, mesg
+ integer :: info
info = -1
- limlev = 6 ! Warning, limlev must cope with grid resolution ...
+ limlev = 2 ! Warning, limlev must cope with grid resolution ...
omega = 1.15_mk ! Relaxation param
- with_w_cycles = .TRUE.
+ with_w_cycles = .FALSE.
dump_info = .TRUE.
ibcvalue(:,:,:,:,:)=0.0_mk
ibcdef(1,:)=bc(:)
@@ -50,12 +62,76 @@ contains
! Anyway, at the time the two parameters above are not used in ppm routine.
!! maxlev : number of levels in multigrid
!!
-
+
call ppm_mg_init(topoid, ppm_param_eq_poisson,ghostsize, ppm_param_smooth_rbsor,dime,ibcdef,&
- ibcvalue,meshid,limlev,with_w_cycles,dump_info,omega,info)
-
+ ibcvalue,meshid,limlev,with_w_cycles,dump_info,omega,info)
+
if(info.ne.0) stop 'Init_multigrid failed ...'
end subroutine init_multigrid
+ subroutine solve_poisson_multigrid(topoid, field, rhs)
+ integer, intent(in) :: topoid
+ real(mk), dimension(:,:,:,:,:), pointer :: field
+ real(mk), dimension(:,:,:,:,:), pointer :: rhs
+
+ integer :: itera, iterf, iter1, iter2, info
+ real(mk) :: Eu
+
+ itera = 3
+ iterf = 3
+ iter1 = 10
+ iter2 = 4
+ print *, 'solve ...', topoid
+
+ ! Bug inside ...
+ call ppm_mg_solv(topoid, field, rhs, dime, itera, iterf, iter1, iter2, Eu, info)
+
+ end subroutine solve_poisson_multigrid
+
+ subroutine init_fftw(fieldin,fieldout, topoid, meshid)
+
+ real(mk), dimension(:,:,:,:,:), pointer :: fieldin, fieldout
+ integer, intent(in) :: topoid, meshid
+
+ ! Parameters for ppm poisson routine
+ integer :: info
+
+ ! Flag to select built-in Green functions (...)
+ integer :: green
+
+ green = ppm_poisson_grn_pois_per ! periodic boundaries
+ info = -1
+ allocate(fftwplanForppm)
+
+ ! Call ppm routine to initialize fftw plan.
+ call ppm_poisson_init(topoid, meshid, fftwplanForppm, fieldin, fieldout, green,info)
+ if(info.NE.0) stop 'PPM Poisson solver init failed.'
+ end subroutine init_fftw
+
+ subroutine solve_poisson_fftw(fieldin, fieldout, topoid, meshid, ghostsize)
+
+ real(mk), dimension(:,:,:,:,:), pointer :: fieldin, fieldout
+ integer, intent(in) :: meshid, topoid
+ integer, dimension(:),pointer :: ghostsize
+
+ integer :: info
+ ! Finite diff scheme used to compute curl
+ integer :: dtype
+
+ info = -1
+ dtype = ppm_poisson_drv_curl_fd2
+
+ ! Solves laplacian(fieldout) = - fieldin
+ call ppm_poisson_solve(topoid, meshid, fftwplanForppm,fieldin,fieldout,ghostsize,info)
+ if(info.NE.0) stop 'PPM Poisson solver failed.'
+
+ info = -1
+ ! Computes fieldout = curl(fieldout) using finite differences, 2nd order (dtype = ppm_poisson_drv_curl_fd2),
+ ! or 4th order (dtype = ppm_poisson_drv_curl_fd4). Last one is untested according to ppm doc ...
+ call ppm_poisson_fd(topoid,meshid,fieldout,fieldout,dtype,info)
+ if(info.NE.0) stop 'PPM Poisson, curl computation failed.'
+
+ end subroutine solve_poisson_fftw
+
end module Solver
diff --git a/HySoP/src/main/Topology.f90 b/HySoP/src/main/Topology.f90
index 0f49bd7975075755ed1f1e31fa52848d66ed7784..68074c446c9982f430efeae9c097b05dedefd3c2 100755
--- a/HySoP/src/main/Topology.f90
+++ b/HySoP/src/main/Topology.f90
@@ -5,7 +5,13 @@ module client_topology
use ppm_module_mktopo
use ppm_module_data, only : ppm_param_assign_internal, ppm_param_decomp_cartesian, ppm_topo
use mpi
- use client_data, only: mk
+ use client_data, only: mk,dime
+
+
+ use ppm_module_numerics_data, only: ppm_param_eq_poisson, ppm_param_smooth_rbsor
+ use ppm_module_mg_init
+
+
implicit none
type(ppm_t_topo), pointer :: topo => null()
@@ -41,9 +47,10 @@ contains
call ppm_mktopo(topoid,meshid, false_xp, false_np, decomposition, assigning, minPos, maxPos, bc, ghostsize, subcost, &
istart, ndata, resolution, info)
if(info.ne.0) stop 'Topology:init initialisation failed'
+
! Get the created topology ...
topo=>ppm_topo(topoid)%t
- mesh=>topo%mesh(meshid)
+ !mesh=>topo%mesh(meshid)
end subroutine init_topo
end module client_topology
diff --git a/HySoP/src/main/main.cxx b/HySoP/src/main/main.cxx
index fcdecbadfb85ba2a39ca952cfe33c1414581d102..ee277ddc157f407015d9962bb5dc15490f80ff6f 100644
--- a/HySoP/src/main/main.cxx
+++ b/HySoP/src/main/main.cxx
@@ -16,13 +16,14 @@ using Parmes::Def::real_t;
extern "C" void createTopoG(int*, int*, int*, double*, double*, int*, double*);
extern "C" void plouhmans();
+extern "C" void NavierStokes3D();
void testPlouhmans()
{
MPI::Init();
assert(MPI::Is_initialized());
- plouhmans();
+ NavierStokes3D();
MPI::Finalize();