Skip to content
Snippets Groups Projects
Commit 3a1ee9e7 authored by Franck Pérignon's avatar Franck Pérignon
Browse files

NS3D with ppm poisson solver (working...?) and partial (draft) penalisation

parent 13807fee
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 471 additions and 74 deletions
HySoP/CMakeLists.txt
+ 1
− 1
View file @ 3a1ee9e7
......@@ -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)
......
HySoP/src/interfaces/ppm/wrap_ppm_topologies.f95
+ 2
− 2
View file @ 3a1ee9e7
......@@ -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()
......
HySoP/src/main/Domain.f90
+ 6
− 3
View file @ 3a1ee9e7
......@@ -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
......
HySoP/src/main/Fields.f90
+ 12
− 39
View file @ 3a1ee9e7
......@@ -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
HySoP/src/main/FieldsComputation.f90 0 → 100755
+ 42
− 0
View file @ 3a1ee9e7
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
HySoP/src/main/NavierStokes3D.f90 0 → 100755
+ 158
− 0
View file @ 3a1ee9e7
!> 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
HySoP/src/main/Penalize.f90 0 → 100755
+ 124
− 0
View file @ 3a1ee9e7
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
HySoP/src/main/Plouhmans.f90
+ 22
− 9
View file @ 3a1ee9e7
......@@ -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
......
HySoP/src/main/Solver.f90
+ 93
− 17
View file @ 3a1ee9e7
......@@ -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
HySoP/src/main/Topology.f90
+ 9
− 2
View file @ 3a1ee9e7
......@@ -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
HySoP/src/main/main.cxx
+ 2
− 1
View file @ 3a1ee9e7
......@@ -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();
......
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