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

Parmes

parent a847df5b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
HySoP/INSTALL
+ 12
− 0
View file @ 99c8f22c
1- create a "build" directory anywhere
2- cmake path-to-sources-of-Parmes -DPPMCore_DIR=/home/perignon/install/ppmcore/share/CMake/ -DPPMNumerics_DIR=/home/perignon/install/ppmnumerics/share/CMake/
3- make -j 8
4- mpirun -np 8 ./parmesRun
HySoP/src/main/Domain.f90
+ 2
− 2
View file @ 99c8f22c
...@@ -40,7 +40,7 @@ contains ...@@ -40,7 +40,7 @@ contains
domain_maxCoords(2) = zFolke domain_maxCoords(2) = zFolke
! 0.3: boundary layer ! 0.3: boundary layer
domain_minCoords(3) = -zFolke - 0.3 domain_minCoords(3) = -zFolke - 0.3
domain_maxCoords(3) = zFolke + 0.3 domain_maxCoords(3) = zFolke + 0.3
! Boundary conditions and ghosts ! Boundary conditions and ghosts
allocate(domain_bc(2*dime), domain_ghostsize(dime), stat = istat) allocate(domain_bc(2*dime), domain_ghostsize(dime), stat = istat)
...@@ -56,7 +56,7 @@ contains ...@@ -56,7 +56,7 @@ contains
allocate(grid_resolution(dime),grid_step(dime),stat=istat) allocate(grid_resolution(dime),grid_step(dime),stat=istat)
if(istat.ne.0) stop 'grid_resolution, allocation failed' if(istat.ne.0) stop 'grid_resolution, allocation failed'
grid_resolution = 65 grid_resolution = 129
grid_step = (domain_maxCoords - domain_minCoords)/(real(grid_resolution,mk)-1.) grid_step = (domain_maxCoords - domain_minCoords)/(real(grid_resolution,mk)-1.)
boundary_layer_min = 0 boundary_layer_min = 0
boundary_layer_max = 0 boundary_layer_max = 0
......
HySoP/src/main/FieldsComputation.f90
+ 16
− 2
View file @ 99c8f22c
...@@ -9,7 +9,7 @@ module FieldsComputation ...@@ -9,7 +9,7 @@ module FieldsComputation
private private
public :: fit_velocity, computeFlowRate public :: fit_velocity, perturb
contains contains
...@@ -55,8 +55,22 @@ contains ...@@ -55,8 +55,22 @@ contains
print *, computed_flowrate, '//',surf print *, computed_flowrate, '//',surf
vel(1,:,:,:,1) = vel(1,:,:,:,1) + constant vel(1,:,:,:,1) = vel(1,:,:,:,1) + constant
end subroutine fit_velocity end subroutine fit_velocity
subroutine perturb(vel, time)
! Velocity, intent(inout)
real(mk), dimension(:,:,:,:,:), pointer :: vel
! current time
real(mk) :: time
real(mk) :: pi
pi = 4.0*atan(1.0_mk)
vel(2,:,:,:,1) = sin(pi*(time-3.0))
end subroutine perturb
subroutine computeFlowRate(vel, stepSize, topo, flowRate) subroutine computeFlowRate(vel, stepSize, topo, flowRate)
......
HySoP/src/main/NavierStokes3D.f90
+ 52
− 35
View file @ 99c8f22c
...@@ -20,7 +20,7 @@ module ppmNavierStokes3D ...@@ -20,7 +20,7 @@ module ppmNavierStokes3D
use penalisation, only : penalisation_init, penalise_velocity, get_drag use penalisation, only : penalisation_init, penalise_velocity, get_drag
use FieldsComputation, only : fit_velocity,computeFlowRate use FieldsComputation, only : fit_velocity,perturb
use vectorcalculus, only: curldf4, computeRHS, norm1, norm2, normInf use vectorcalculus, only: curldf4, computeRHS, norm1, norm2, normInf
...@@ -55,7 +55,7 @@ contains ...@@ -55,7 +55,7 @@ contains
! debug mode ! debug mode
integer :: debug integer :: debug
! error status ! error status
integer :: info integer :: info,rank,sizeZ
!====================== !======================
! Init ppm ! Init ppm
...@@ -82,6 +82,7 @@ contains ...@@ -82,6 +82,7 @@ contains
! Fields allocation ! Fields allocation
!====================== !======================
call init_fields(domain_ghostsize, topo) call init_fields(domain_ghostsize, topo)
!====================== !======================
! Init Physics ! Init Physics
!====================== !======================
...@@ -94,12 +95,30 @@ contains ...@@ -94,12 +95,30 @@ contains
velocity = 0.0_mk velocity = 0.0_mk
call init_vorticity(vorticity,topo,grid_step,boundary_layer_min,boundary_layer_max) call init_vorticity(vorticity,topo,grid_step,boundary_layer_min,boundary_layer_max)
call penalisation_init(topo,grid_step,boundary_layer_min,boundary_layer_max) call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info)
!!$ vorticity = 0.0_mk
!!$ if(rank==0) then
!!$ sizeZ = size(vorticity,4)
!!$ print *, 'shape : ', shape(vorticity), sizeZ
!!$ vorticity(2,4,5,5,:) = 12.0
!!$ else
!!$ vorticity(2,4,5,1,:) = 12.0
!!$ end if
!!$ call MPI_BARRIER(MPI_COMM_WORLD,info)
!====================== call ppm_map_field_ghost_get(topo%ID,topo%mesh(1)%ID,domain_ghostsize,info)
! Init Particles call ppm_map_field_push(topo%ID,topo%mesh(1)%ID,vorticity, 3, info)
!====================== call ppm_map_field_send(info)
call ppm_map_field_pop(topo%ID,topo%mesh(1)%ID,vorticity,3,domain_ghostsize,info)
!============================================
! Init penalization : set chi functions
!============================================
call penalisation_init(topo,grid_step,boundary_layer_min,boundary_layer_max)
!============================================
! Build a first set of particles
!============================================
call init_particles(vorticity,topo%ID,topo%mesh(1)%ID,velocity) call init_particles(vorticity,topo%ID,topo%mesh(1)%ID,velocity)
print *, "end of parmes:ppm:init_client" print *, "end of parmes:ppm:init_client"
...@@ -120,12 +139,14 @@ contains ...@@ -120,12 +139,14 @@ contains
Re = 133.4 Re = 133.4
nu = 1./Re nu = 1./Re
initial_time = 0.0 initial_time = 0.0
final_time = 10.0 final_time = 60.0
time_step = 0.05 time_step = 0.04
current_time = initial_time current_time = initial_time
iter = 1 iter = 1
max_iter= 10 max_iter= 10000
! call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info)
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info)
! call MPI_COMM_SIZE(MPI_COMM_WORLD,nbProcs,info) ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nbProcs,info)
! Initialisation stuff ! Initialisation stuff
call init_client() call init_client()
...@@ -136,17 +157,9 @@ contains ...@@ -136,17 +157,9 @@ contains
! Init multigrid works in ppm (or seems to ...) but solve fails. ! 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_multigrid(topo%ID, topo%mesh(1)%ID, domain_ghostsize, domain_bc)
call MPI_BARRIER(MPI_COMM_WORLD,info)
! Analytical solutions ! Analytical solutions
!call rhs_analytic(topo, grid_step, stream_function,vel_ex,vorticity,nu) !call rhs_analytic(topo, grid_step, stream_function,vel_ex,vorticity,nu)
! Init ghosts ...
call ppm_map_field_ghost_get(topo%ID,topo%mesh(1)%ID,domain_ghostsize,info)
call ppm_map_field_push(topo%ID,topo%mesh(1)%ID,vorticity, 3, info)
call ppm_map_field_send(info)
call ppm_map_field_pop(topo%ID,topo%mesh(1)%ID,vorticity,3,domain_ghostsize,info)
derive = ppm_poisson_drv_curl_fd4 derive = ppm_poisson_drv_curl_fd4
call init_poisson_solver(vorticity,velocity,topo%ID, topo%mesh(1)%ID,derive) call init_poisson_solver(vorticity,velocity,topo%ID, topo%mesh(1)%ID,derive)
...@@ -155,13 +168,15 @@ contains ...@@ -155,13 +168,15 @@ contains
print *, '[',rank,'] pas espace/temps ', grid_step, time_step print *, '[',rank,'] pas espace/temps ', grid_step, time_step
print *, '[',rank,'] mesh ', shape(topo%mesh(1)%nnodes), topo%mesh(1)%nnodes print *, '[',rank,'] mesh ', shape(topo%mesh(1)%nnodes), topo%mesh(1)%nnodes
allocate(drag(max_iter)) !if(rank==0) then
allocate(drag(max_iter))
open(10,file='drag') open(10,file='drag')
!end if
call MPI_BARRIER(MPI_COMM_WORLD,info)
print *, 'start iterations ...'
! Time loop ! Time loop
do while(current_time < final_time .and. iter <= max_iter) do while(current_time < final_time .and. iter <= max_iter)
print *, '================ Start iter ',iter,' ================' ! print *, '================ Start iter ',iter,' ================'
! Compute velocity from vorticity ! Compute velocity from vorticity
! Two steps: ! Two steps:
...@@ -169,22 +184,24 @@ contains ...@@ -169,22 +184,24 @@ contains
! to be done in ppm routine; indeed we do not have to deal with stream_function. ! to be done in ppm routine; indeed we do not have to deal with stream_function.
! - update velocity to fit with required boundary values ! - update velocity to fit with required boundary values
call solve_poisson(vorticity, velocity, topo%ID, topo%mesh(1)%ID,domain_ghostsize) call solve_poisson(vorticity, velocity, topo%ID, topo%mesh(1)%ID,domain_ghostsize)
call MPI_BARRIER(MPI_COMM_WORLD,info) call fit_velocity(velocity, grid_step,lengths,topo,boundary_layer_min,boundary_layer_max)
!call fit_velocity(velocity, grid_step,lengths,topo,boundary_layer_min,boundary_layer_max) if( current_time>3.0 .and. current_time<4.0) then
call perturb(velocity, current_time)
end if
! Penalize velocity ! Penalize velocity
call penalise_velocity(velocity,topo,time_step) call penalise_velocity(velocity,topo,time_step)
call ppm_map_field_ghost_get(topo%ID,topo%mesh(1)%ID,domain_ghostsize,info) call ppm_map_field_ghost_get(topo%ID,topo%mesh(1)%ID,domain_ghostsize,info)
call ppm_map_field_push(topo%ID,topo%mesh(1)%ID,velocity, 3, info) call ppm_map_field_push(topo%ID,topo%mesh(1)%ID,velocity, 3, info)
call ppm_map_field_send(info) call ppm_map_field_send(info)
call ppm_map_field_pop(topo%ID,topo%mesh(1)%ID,velocity,3,domain_ghostsize,info) call ppm_map_field_pop(topo%ID,topo%mesh(1)%ID,velocity,3,domain_ghostsize,info)
call get_drag(drag(iter))
!print *, 'drag...', drag(iter,rank)
write(10,'(e14.5,e14.5)') current_time, drag(iter)
!!$ ! Compute the new "penalized" vorticity call get_drag(drag(iter))
if(rank == 0) then
write(10,'(e14.5,e14.5)') current_time, drag(iter)
end if
! Compute the new "penalized" vorticity
call curlDF4(velocity, vorticity,topo, grid_step) call curlDF4(velocity, vorticity,topo, grid_step)
! Computation of rhs of vorticity equation (on the grid) ... ! Computation of rhs of vorticity equation (on the grid) ...
call computeRHS(velocity, vorticity, rhs, topo, grid_step, nu) call computeRHS(velocity, vorticity, rhs, topo, grid_step, nu)
call MPI_BARRIER(MPI_COMM_WORLD,info) call MPI_BARRIER(MPI_COMM_WORLD,info)
...@@ -206,12 +223,13 @@ contains ...@@ -206,12 +223,13 @@ contains
! Remesh ! Remesh
call remesh(topo%ID, topo%mesh(1)%ID,domain_ghostsize,vorticity) call remesh(topo%ID, topo%mesh(1)%ID,domain_ghostsize,vorticity)
! call printToVTK("post",iter, current_time,velocity,vorticity,topo,grid_step) if(mod(iter, 10).eq.0) then
call printToVTK("post",iter, current_time,velocity,vorticity,topo,grid_step)
end if
current_time = current_time + time_step current_time = current_time + time_step
!current_time = final_time !current_time = final_time
iter = iter+1 iter = iter+1
print *, 'end of iteration ', iter-1 !print *, 'end of iteration ', iter-1
call MPI_BARRIER(MPI_COMM_WORLD,info) call MPI_BARRIER(MPI_COMM_WORLD,info)
enddo enddo
...@@ -219,7 +237,6 @@ contains ...@@ -219,7 +237,6 @@ contains
close(10) close(10)
call ppm_finalize(info) call ppm_finalize(info)
print *, 'end ppm simulation' print *, 'end ppm simulation'
deallocate(drag) deallocate(drag)
end subroutine main_client end subroutine main_client
......
HySoP/src/main/Particles.f90
+ 107
− 44
View file @ 99c8f22c
...@@ -2,9 +2,9 @@ ...@@ -2,9 +2,9 @@
module Particles module Particles
use ppm_module_rmsh, only : ppm_rmsh_create_part,ppm_interp_m2p, ppm_interp_p2m !, ppm_rmsh_remesh use ppm_module_rmsh, only : ppm_rmsh_create_part,ppm_interp_m2p, ppm_interp_p2m !, ppm_rmsh_remesh
use ppm_module_impose_part_bc
use client_data, only: mk, dime use client_data, only: mk, dime
! use mpi use mpi
use ppm_module_data, only : ppm_param_rmsh_kernel_mp4 use ppm_module_data, only : ppm_param_rmsh_kernel_mp4
use ppm_module_map_part use ppm_module_map_part
use ppm_module_util_dbg use ppm_module_util_dbg
...@@ -47,7 +47,7 @@ contains ...@@ -47,7 +47,7 @@ contains
! topo and mesh ids ! topo and mesh ids
integer, intent(in) :: topoid integer, intent(in) :: topoid
integer, intent(in) :: meshid integer, intent(in) :: meshid
integer :: info, newNpart integer :: info, newNpart,rank,i
! velocity on grid ! velocity on grid
real(mk), dimension(:,:,:,:,:), pointer :: vel real(mk), dimension(:,:,:,:,:), pointer :: vel
...@@ -62,31 +62,66 @@ contains ...@@ -62,31 +62,66 @@ contains
npart = 0 npart = 0
last_npart = npart last_npart = npart
first_call = .TRUE. first_call = .TRUE.
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info)
! --------> The following call will allocate memory for xp, omp and velo_p, right? ! --------> The following call will allocate memory for xp, omp and velo_p, right?
! --------> And deallocation will be handled by ppm also? ! --------> And deallocation will be handled by ppm also?
call ppm_rmsh_create_part(topoid,meshid,xp,npart,omp,dime,field_on_grid,cutoff,info,.TRUE.,& call ppm_rmsh_create_part(topoid,meshid,xp,npart,omp,dime,field_on_grid,cutoff,info,.TRUE.,&
field_wp=vel,wp=velo_p,lda2=dime) field_wp=vel,wp=velo_p,lda2=dime)
print *, '[',rank,'] initialisation with ', npart,' particles'
!> --------> Other variables carried by particles but not allocated during ppm_rmsh_create call !> --------> Other variables carried by particles but not allocated during ppm_rmsh_create call
! --------> is there a better way to do this, with map_push or other ppm routines? ! --------> is there a better way to do this, with map_push or other ppm routines?
allocate(rhsp(dime,npart),xp0(dime,npart)) allocate(rhsp(dime,npart),xp0(dime,npart))
!call ppm_dbg_print_d(topoid,0.0_mk,0,3,info,xp,npart)
! --------> Is the following (map_part_global) required? I would say yes but I'm not sure ...
newNpart = 0 ! print *, '[',rank,'] shapes ...', shape(xp),'//', shape(omp), '//', shape(velo_p),'//',shape(rhsp),'//', shape(xp0)
call ppm_map_part_global(topoid,xp,npart,info) ! positions
call ppm_map_part_push(velo_p,dime,npart,info) ! velocity
call ppm_map_part_push(rhsp,dime,npart,info) ! rhs
call ppm_map_part_push(omp,dime,npart,info) ! vorticity
call ppm_map_part_push(xp0,dime,npart,info)
call ppm_map_part_send(npart,newNpart,info) ! send
call ppm_map_part_pop(xp0,dime,npart,newNpart,info)
call ppm_map_part_pop(omp,dime,npart,newNpart,info)
call ppm_map_part_pop(velo_p,dime,npart,newNpart,info)
call ppm_map_part_pop(rhsp,dime,npart,newNpart,info)
call ppm_map_part_pop(xp,dime,npart,newNpart,info)
npart = newNpart
!!$ ! --------> Is the following (map_part_global) required? I would say no but I'm not sure ...
!!$ newNpart = 0
!!$ call ppm_impose_part_bc(topoid,xp,npart,info)
!!$ call ppm_map_part_global(topoid,xp,npart,info) ! positions
!!$ call ppm_map_part_push(velo_p,dime,npart,info) ! velocity
!!$ call ppm_map_part_push(rhsp,dime,npart,info) ! rhs
!!$ call ppm_map_part_push(omp,dime,npart,info) ! vorticity
!!$ call ppm_map_part_push(xp0,dime,npart,info)
!!$ call ppm_map_part_send(npart,newNpart,info) ! send
!!$ call ppm_map_part_pop(xp0,dime,npart,newNpart,info)
!!$ call ppm_map_part_pop(omp,dime,npart,newNpart,info)
!!$ call ppm_map_part_pop(rhsp,dime,npart,newNpart,info)
!!$ call ppm_map_part_pop(velo_p,dime,npart,newNpart,info)
!!$ call ppm_map_part_pop(xp,dime,npart,newNpart,info)
!!$
!!$ npart = newNpart
! print *, '[',rank,'] post init with ', npart,' particles'
! print *, '[',rank,'] post init xp/omp ', xp(3,:),velo_p(:,:),omp
!!$ if(rank==1) then
!!$ do i =1,size(xp,2)
!!$ xp(3,i) = xp(3,i) - 1.0_mk
!!$ end do
!!$ else
!!$ xp(3,:) = xp(3,:) +1.5
!!$ end if
!!$ !call ppm_impose_part_bc(topoid,xp,npart,info)
!!$ call ppm_map_part_partial(topoid,xp,npart,info) ! positions
!!$ call ppm_map_part_push(velo_p,dime,npart,info) ! velocity
!!$ !call ppm_map_part_push(rhsp,dime,npart,info) ! rhs
!!$ !call ppm_map_part_push(omp,dime,npart,info) ! vorticity
!!$ !call ppm_map_part_push(xp0,dime,npart,info)
!!$ call ppm_map_part_send(npart,newNpart,info) ! send
!!$ !call ppm_map_part_pop(xp0,dime,npart,newNpart,info)
!!$ !call ppm_map_part_pop(omp,dime,npart,newNpart,info)
!!$ !call ppm_map_part_pop(rhsp,dime,npart,newNpart,info)
!!$ call ppm_map_part_pop(velo_p,dime,npart,newNpart,info)
!!$ call ppm_map_part_pop(xp,dime,npart,newNpart,info)
!!$ !call ppm_impose_part_bc(topoid,xp,npart,info)
!!$ npart = newNpart
!!$ print *, '[',rank,'] post 2 init with ', npart,' particles'
!!$ print *, '[',rank,'] post 2 init xp/omp ', xp(3,:)!,velo_p(:,:)
!!$
!!$
!!$
end subroutine init_particles end subroutine init_particles
!> Create particles distribution !> Create particles distribution
...@@ -98,7 +133,7 @@ contains ...@@ -98,7 +133,7 @@ contains
! topo and mesh ids ! topo and mesh ids
integer, intent(in) :: topoid integer, intent(in) :: topoid
integer, intent(in) :: meshid integer, intent(in) :: meshid
integer :: info, newNpart integer :: info, newNpart,rank
integer, dimension(:), pointer :: ghostsize integer, dimension(:), pointer :: ghostsize
! velocity on grid ! velocity on grid
real(mk), dimension(:,:,:,:,:), pointer :: vel real(mk), dimension(:,:,:,:,:), pointer :: vel
...@@ -107,27 +142,33 @@ contains ...@@ -107,27 +142,33 @@ contains
info = 0 info = 0
!call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info) call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info)
! Call ppm func to create the particles distribution ! Call ppm func to create the particles distribution
last_npart = npart last_npart = npart
call ppm_rmsh_create_part(topoid,meshid,xp,npart,omp,dime,field_on_grid,cutoff,info,resetpos,& call ppm_rmsh_create_part(topoid,meshid,xp,npart,omp,dime,field_on_grid,cutoff,info,resetpos,&
field_wp=vel,wp=velo_p,lda2=dime) field_wp=vel,wp=velo_p,lda2=dime)
!if(associated(xp)) print *, "xp shape:", rank, " ", shape(xp), npart !if(associated(xp)) print *, "xp shape:", rank, " ", shape(xp), npart
if(last_npart.lt.npart) then
newNpart = 0 deallocate(rhsp,xp0)
call ppm_map_part_partial(topoid,xp,npart,info) ! positions
call ppm_map_part_push(velo_p,dime,npart,info) ! velocity
call ppm_map_part_send(npart,newNpart,info) ! send
call ppm_map_part_pop(velo_p,dime,npart,newNpart,info)
call ppm_map_part_pop(xp,dime,npart,newNpart,info)
npart = newNpart
! --------> rhsp and xp0 size was computed before rmsh_create and part_pop => must be updated
! --------> deallocate/allocate may be done by some ppm routine?
if (last_npart.lt.npart) then
deallocate(rhsp,xp0) ! Todo : use ppm_alloc stuff
allocate(rhsp(dime,npart),xp0(dime,npart)) allocate(rhsp(dime,npart),xp0(dime,npart))
end if end if
! newNpart = 0
!!$ call ppm_impose_part_bc(topoid,xp,npart,info)
!!$ call ppm_map_part_global(topoid,xp,npart,info) ! positions
!!$ call ppm_map_part_push(velo_p,dime,npart,info) ! velocity
!!$ call ppm_map_part_push(rhsp,dime,npart,info) ! rhs
!!$ call ppm_map_part_push(omp,dime,npart,info) ! vorticity
!!$ call ppm_map_part_push(xp0,dime,npart,info)
!!$ call ppm_map_part_send(npart,newNpart,info) ! send
!!$ call ppm_map_part_pop(xp0,dime,npart,newNpart,info)
!!$ call ppm_map_part_pop(omp,dime,npart,newNpart,info)
!!$ call ppm_map_part_pop(rhsp,dime,npart,newNpart,info)
!!$ call ppm_map_part_pop(velo_p,dime,npart,newNpart,info)
!!$ call ppm_map_part_pop(xp,dime,npart,newNpart,info)
!print *, 'UPDATE ', npart, newNpart, shape(omp) , '\\', shape(velo_p)
! npart = newNpart
! --------> Ok, from here all vars carried by particles must have the following shape : dime,npart ! --------> Ok, from here all vars carried by particles must have the following shape : dime,npart
...@@ -136,7 +177,7 @@ contains ...@@ -136,7 +177,7 @@ contains
!> --------> mesh to particles for rhs of vorticity eq. This routine does not resize rhsp, it must be of the right size before the call? !> --------> mesh to particles for rhs of vorticity eq. This routine does not resize rhsp, it must be of the right size before the call?
call ppm_interp_m2p(topoid, meshid, xp, npart, rhsp, dime, kernel, ghostsize, rhs,info) call ppm_interp_m2p(topoid, meshid, xp, npart, rhsp, dime, kernel, ghostsize, rhs,info)
print *, 'End of update particles' ! print *, 'End of update particles'
end subroutine update_particles end subroutine update_particles
...@@ -156,8 +197,8 @@ contains ...@@ -156,8 +197,8 @@ contains
real(mk), dimension(:,:,:,:,:), pointer :: rhs real(mk), dimension(:,:,:,:,:), pointer :: rhs
!> local loop indices !> local loop indices
integer :: i, newNpart integer :: i, newNpart,rank
! call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info) call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info)
do i=1,npart do i=1,npart
xp0(1,i) = xp(1,i) xp0(1,i) = xp(1,i)
...@@ -170,7 +211,9 @@ contains ...@@ -170,7 +211,9 @@ contains
end do end do
newNpart = 0 newNpart = 0
! Particles positions have changed ... we must map between domains
! call ppm_impose_part_bc(topoid,xp,npart,info) ! Useless according to doc, required according to Omar ...
call ppm_map_part_partial(topoid,xp,npart,info) ! positions call ppm_map_part_partial(topoid,xp,npart,info) ! positions
call ppm_map_part_push(velo_p,dime,npart,info) ! velocity call ppm_map_part_push(velo_p,dime,npart,info) ! velocity
call ppm_map_part_push(rhsp,dime,npart,info) ! rhs call ppm_map_part_push(rhsp,dime,npart,info) ! rhs
...@@ -179,10 +222,17 @@ contains ...@@ -179,10 +222,17 @@ contains
call ppm_map_part_send(npart,newNpart,info) ! send call ppm_map_part_send(npart,newNpart,info) ! send
call ppm_map_part_pop(xp0,dime,npart,newNpart,info) call ppm_map_part_pop(xp0,dime,npart,newNpart,info)
call ppm_map_part_pop(omp,dime,npart,newNpart,info) call ppm_map_part_pop(omp,dime,npart,newNpart,info)
call ppm_map_part_pop(velo_p,dime,npart,newNpart,info) call ppm_map_part_pop(rhsp,dime,npart,newNpart,info)
call ppm_map_part_pop(rhsp,dime,npart,newNpart,info) call ppm_map_part_pop(velo_p,dime,npart,newNpart,info)
call ppm_map_part_pop(xp,dime,npart,newNpart,info) call ppm_map_part_pop(xp,dime,npart,newNpart,info)
! print *, 'PUSH 0 ', npart, newNpart, shape(omp) , '\\', shape(velo_p)
! Update sizes ...
! if (newNpart.lt.npart) then
! deallocate(rhsp,xp0) ! Todo : use ppm_alloc stuff
! allocate(rhsp(dime,newNpart),xp0(dime,newNpart))
! end if
npart = newNpart npart = newNpart
! --------> No need to call any realloc or resize process for xp0, rhsp ... that was done during map_part_pop, right? ! --------> No need to call any realloc or resize process for xp0, rhsp ... that was done during map_part_pop, right?
...@@ -203,11 +253,24 @@ contains ...@@ -203,11 +253,24 @@ contains
end do end do
newNpart = 0 newNpart = 0
!call ppm_impose_part_bc(topoid,xp,npart,info)
call ppm_map_part_partial(topoid,xp,npart,info) ! positions call ppm_map_part_partial(topoid,xp,npart,info) ! positions
call ppm_map_part_push(omp,dime,npart,info) ! strengths call ppm_map_part_push(velo_p,dime,npart,info) ! velocity
call ppm_map_part_push(rhsp,dime,npart,info) ! rhs
call ppm_map_part_push(omp,dime,npart,info) ! vorticity
call ppm_map_part_push(xp0,dime,npart,info)
call ppm_map_part_send(npart,newNpart,info) ! send call ppm_map_part_send(npart,newNpart,info) ! send
call ppm_map_part_pop(omp,dime,npart,newNpart,info) ! strengths call ppm_map_part_pop(xp0,dime,npart,newNpart,info)
call ppm_map_part_pop(xp,dime,npart,newNpart,info) ! positions call ppm_map_part_pop(omp,dime,npart,newNpart,info)
call ppm_map_part_pop(rhsp,dime,npart,newNpart,info)
call ppm_map_part_pop(velo_p,dime,npart,newNpart,info)
call ppm_map_part_pop(xp,dime,npart,newNpart,info)
! print *, 'PUSH 1 ', npart, newNpart, shape(omp) , '\\', shape(velo_p)
! if (newNpart.lt.npart) then
! deallocate(rhsp,xp0) ! Todo : use ppm_alloc stuff
! allocate(rhsp(dime,newNpart),xp0(dime,newNpart))
! end if
npart = newNpart npart = newNpart
end subroutine push_particles end subroutine push_particles
......
HySoP/src/main/Penalize.f90
+ 5
− 3
View file @ 99c8f22c
...@@ -94,13 +94,14 @@ contains ...@@ -94,13 +94,14 @@ contains
chi_sphere(i,j,k,isub) = lambda chi_sphere(i,j,k,isub) = lambda
end if end if
if( (z> boundary_layer_max(3)).or.(z< boundary_layer_min(3))) then if( (z> boundary_layer_max(3)).or.(z< boundary_layer_min(3))) then
chi(i,j,k,isub) = lambda chi(i,j,k,isub) = lambda
end if end if
chi(i,j,k,isub) = chi(i,j,k,isub) + chi_sphere(i,j,k,isub)
end do end do
end do end do
end do end do
end do end do
chi = chi + chi_sphere
nullify(isublist, ndata) nullify(isublist, ndata)
local_drag = 0.0 local_drag = 0.0
...@@ -149,11 +150,12 @@ subroutine get_drag(drag) ...@@ -149,11 +150,12 @@ subroutine get_drag(drag)
real(mk), intent(out) :: drag real(mk), intent(out) :: drag
integer :: info integer :: info
real(mk) :: coef,uinf real(mk) :: coef,uinf
! call MPI_BARRIER(MPI_COMM_WORLD,info)
uinf = 1.0 uinf = 1.0
coef = 2./(uinf**2*2.*sphere_radius) coef = 2./(uinf**2*2.*sphere_radius)
local_drag = coef*local_drag local_drag = coef*local_drag
call MPI_AllReduce(local_drag, drag, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD,info) call MPI_AllReduce(local_drag, drag, 1, MPI_DOUBLE_PRECISION, MPI_SUM,MPI_COMM_WORLD,info)
end subroutine get_drag end subroutine get_drag
......
HySoP/src/main/Topology.f90
+ 3
− 3
View file @ 99c8f22c
...@@ -4,7 +4,7 @@ module client_topology ...@@ -4,7 +4,7 @@ module client_topology
use ppm_module_typedef, only : ppm_t_topo use ppm_module_typedef, only : ppm_t_topo
use ppm_module_mktopo use ppm_module_mktopo
use ppm_module_data, only : ppm_param_assign_internal, ppm_param_decomp_cartesian, ppm_topo use ppm_module_data, only : ppm_param_assign_internal, ppm_param_decomp_cartesian, ppm_topo
! use mpi use mpi
use client_data, only: mk,dime use client_data, only: mk,dime
...@@ -48,7 +48,7 @@ contains ...@@ -48,7 +48,7 @@ contains
istart, ndata, resolution, info) istart, ndata, resolution, info)
if(info.ne.0) stop 'Topology:init initialisation failed' if(info.ne.0) stop 'Topology:init initialisation failed'
! call MPI_COMM_RANK(MPI_COMM_WORLD,rank,info) call MPI_COMM_RANK(MPI_COMM_WORLD,rank,info)
! print *, '[',rank,'], topo :', shape(ndata), '//', ndata ! print *, '[',rank,'], topo :', shape(ndata), '//', ndata
! Get the created topology ... ! Get the created topology ...
topo=>ppm_topo(topoid)%t topo=>ppm_topo(topoid)%t
...@@ -61,7 +61,7 @@ contains ...@@ -61,7 +61,7 @@ contains
! sub2proc[i] = global number of the processus on which subdomain i is attached ! sub2proc[i] = global number of the processus on which subdomain i is attached
! !
!!$ !!$
print *, '[', rank,'], subs : ', topo%nsublist, '//', shape(topo%isublist),'//', topo%isublist print *, '[', rank,'], subs : ', topo%nsublist, '//', shape(topo%isublist),'//', topo%isublist(1)
!!$ print *, '======================= [', rank,'] ', shape(topo%sub2proc) !!$ print *, '======================= [', rank,'] ', shape(topo%sub2proc)
!!$ print *, '======================= [', rank,'] ', topo%sub2proc !!$ print *, '======================= [', rank,'] ', topo%sub2proc
!!$ !!$
......
HySoP/src/main/io_vtk.f90
+ 3
− 5
View file @ 99c8f22c
...@@ -3,11 +3,9 @@ module io_vtk ...@@ -3,11 +3,9 @@ module io_vtk
use client_data, only:mk use client_data, only:mk
use client_topology, only: ppm_t_topo use client_topology, only: ppm_t_topo
!use mpi use mpi
implicit none implicit none
contains contains
!> print velocity and vorticity to vtk file !> print velocity and vorticity to vtk file
...@@ -30,7 +28,7 @@ contains ...@@ -30,7 +28,7 @@ contains
integer :: isub,isubl, rank integer :: isub,isubl, rank
! local resolution ! local resolution
integer :: nx,ny,nz integer :: nx,ny,nz
! integer :: info integer :: info
!local filename (depends on subproc number) !local filename (depends on subproc number)
character(len=30) :: localname,buffer1, buffer2,localname2 character(len=30) :: localname,buffer1, buffer2,localname2
...@@ -45,7 +43,7 @@ contains ...@@ -45,7 +43,7 @@ contains
!call res_vtk_vit(trim(name_vit)) !call res_vtk_vit(trim(name_vit))
!call res_vtk_omega(trim(name_omg)) !call res_vtk_omega(trim(name_omg))
!call MPI_COMM_RANK(MPI_COMM_WORLD,rank,info) call MPI_COMM_RANK(MPI_COMM_WORLD,rank,info)
do isub=1,topo%nsublist do isub=1,topo%nsublist
isubl=topo%isublist(isub) isubl=topo%isublist(isub)
......
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