From bd55acc12a82c5bca9f21ddded2e5c1b8fbf5de0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Wed, 13 Jul 2011 15:07:31 +0000
Subject: [PATCH] Add lift computation
---
HySoP/INSTALL | 7 +++
HySoP/src/main/Domain.f90 | 3 +-
HySoP/src/main/FieldsComputation.f90 | 19 +++----
HySoP/src/main/NavierStokes3D.f90 | 58 ++++++++++++--------
HySoP/src/main/Particles.f90 | 31 ++---------
HySoP/src/main/Penalize.f90 | 26 ++++++---
HySoP/src/main/TestFunctions.f90 | 5 +-
HySoP/src/main/client_data.f90 | 82 +---------------------------
8 files changed, 79 insertions(+), 152 deletions(-)
diff --git a/HySoP/INSTALL b/HySoP/INSTALL
index ede41e12c..6db0240a3 100644
--- a/HySoP/INSTALL
+++ b/HySoP/INSTALL
@@ -9,4 +9,11 @@
4- mpirun -np 8 ./parmesRun
+Options :
+
+debug -> cmake -DCMAKE_BUILD_TYPE=Debug
+
+verbose-> make VERBOSE=1
+
+
diff --git a/HySoP/src/main/Domain.f90 b/HySoP/src/main/Domain.f90
index 2487885d4..8dcae7927 100755
--- a/HySoP/src/main/Domain.f90
+++ b/HySoP/src/main/Domain.f90
@@ -33,7 +33,6 @@ contains
allocate(boundary_layer_min(dime), boundary_layer_max(dime), stat = istat)
if(istat.ne.0) stop 'Geometry, allocation failed'
-
domain_minCoords(1) = 0.0
domain_maxCoords(1) = 2.*zFolke
domain_minCoords(2) = -zFolke
@@ -56,7 +55,7 @@ contains
allocate(grid_resolution(dime),grid_step(dime),stat=istat)
if(istat.ne.0) stop 'grid_resolution, allocation failed'
- grid_resolution = 129
+ grid_resolution = 65
grid_step = (domain_maxCoords - domain_minCoords)/(real(grid_resolution,mk)-1.)
boundary_layer_min = 0
boundary_layer_max = 0
diff --git a/HySoP/src/main/FieldsComputation.f90 b/HySoP/src/main/FieldsComputation.f90
index 1803b5bbd..1802e0ae7 100755
--- a/HySoP/src/main/FieldsComputation.f90
+++ b/HySoP/src/main/FieldsComputation.f90
@@ -1,6 +1,6 @@
module FieldsComputation
- use client_data, only: mk
+ use client_data, only: mk,pi
use client_topology, only: ppm_t_topo
use mpi
!use Solver, only: solve_poisson
@@ -44,16 +44,15 @@ contains
! The flow is supposed to be along direction 1
! physical domain surface
- surf = (boundary_layer_max(3)-boundary_layer_min(3))*lengths(2)
+ !surf = (boundary_layer_max(3)-boundary_layer_min(3))*lengths(2)
! Value of the flow rate we want to fix
- theoretical_flowrate = 1.
- call computeFlowRate(vel, stepSize, topo,computed_flowrate)
+ !theoretical_flowrate = 1.
+ !call computeFlowRate(vel, stepSize, topo,computed_flowrate)
! Correction is added to velocity ...
constant = (boundary_layer_max(3)-boundary_layer_min(3))/lengths(3) !(theoretical_flowrate - computed_flowrate)/surf
-! print *, 'const', constant
-
- ! print *, computed_flowrate, '//',surf
-
+ ! print *, 'const', constant
+ ! print *, computed_flowrate, '//',surf
+
vel(1,:,:,:,1) = vel(1,:,:,:,1) + constant
end subroutine fit_velocity
@@ -63,10 +62,6 @@ contains
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))
diff --git a/HySoP/src/main/NavierStokes3D.f90 b/HySoP/src/main/NavierStokes3D.f90
index 2c872d60f..2d2abae41 100755
--- a/HySoP/src/main/NavierStokes3D.f90
+++ b/HySoP/src/main/NavierStokes3D.f90
@@ -8,7 +8,7 @@ module ppmNavierStokes3D
use ppm_module_map_field
! use client_io
- use client_data, only: mk, dime
+ use client_data, only: mk, dime,rank
! Physical domain and grid
use Domain
! Fields on the grid
@@ -18,7 +18,7 @@ module ppmNavierStokes3D
! Multigrid solver
use Solver, only : init_poisson_solver, solve_poisson, ppm_poisson_drv_none, ppm_poisson_drv_curl_fd2, ppm_poisson_drv_curl_fd4
- use penalisation, only : penalisation_init, penalise_velocity, get_drag
+ use penalisation, only : penalisation_init, penalise_velocity, get_dragAndLift
use FieldsComputation, only : fit_velocity,perturb
@@ -55,7 +55,7 @@ contains
! debug mode
integer :: debug
! error status
- integer :: info,rank,sizeZ
+ integer :: info,sizeZ
!======================
! Init ppm
@@ -128,25 +128,28 @@ contains
subroutine main_client() bind(c,name='NavierStokes3D')
real(mk) :: initial_time, final_time, current_time, time_step
- integer :: rank,iter, derive,max_iter!,nbProcs
+ integer :: iter, derive,max_iter!,nbProcs
logical :: resetpos
real(mk) :: nu,Re
! real(mk), dimension(:,:), pointer :: err1, err2, errinf
real(mk), dimension(:), pointer :: drag =>NULL()
+ real(mk), dimension(:), pointer :: lifty =>NULL()
+ real(mk), dimension(:), pointer :: liftz =>NULL()
- Re = 133.4
+ real(mk) :: fit_velo
+
+ Re = 210.6 !133.4
nu = 1./Re
initial_time = 0.0
- final_time = 60.0
- time_step = 0.04
+ final_time = 100.0
+ time_step = 0.2
current_time = initial_time
iter = 1
- max_iter= 10000
-
+ max_iter= int((final_time -initial_time)/time_step) + 1
+ print *, 'aaaaaaaaa', max_iter
- call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info)
! call MPI_COMM_SIZE(MPI_COMM_WORLD,nbProcs,info)
! Initialisation stuff
call init_client()
@@ -168,15 +171,18 @@ contains
print *, '[',rank,'] pas espace/temps ', grid_step, time_step
print *, '[',rank,'] mesh ', shape(topo%mesh(1)%nnodes), topo%mesh(1)%nnodes
- !if(rank==0) then
- allocate(drag(max_iter))
- open(10,file='drag')
- !end if
+ allocate(drag(max_iter),lifty(max_iter),liftz(max_iter))
+ if(rank==0) open(10,file='drag')
+
+ fit_velo = (boundary_layer_max(3)-boundary_layer_min(3))/lengths(3)
+
+ ! Synchronise all procs after initialization process ...
call MPI_BARRIER(MPI_COMM_WORLD,info)
+
print *, 'start iterations ...'
! Time loop
do while(current_time < final_time .and. iter <= max_iter)
-! print *, '================ Start iter ',iter,' ================'
+ ! print *, '================ Start iter ',iter,' ================'
! Compute velocity from vorticity
! Two steps:
@@ -184,27 +190,35 @@ contains
! to be done in ppm routine; indeed we do not have to deal with 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_step,lengths,topo,boundary_layer_min,boundary_layer_max)
+ !call fit_velocity(velocity, grid_step,lengths,topo,boundary_layer_min,boundary_layer_max)
+
+ velocity(1,:,:,:,1) = velocity(1,:,:,:,1) + fit_velo
+
+ ! Perturbation step, only required for high Reynolds number
+
if( current_time>3.0 .and. current_time<4.0) then
call perturb(velocity, current_time)
end if
- ! Penalize velocity
+
+ ! Penalize velocity and map result to all procs
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_push(topo%ID,topo%mesh(1)%ID,velocity, 3, info)
call ppm_map_field_send(info)
call ppm_map_field_pop(topo%ID,topo%mesh(1)%ID,velocity,3,domain_ghostsize,info)
- call get_drag(drag(iter))
+ ! The drag ...
+ call get_dragAndLift(drag(iter),lifty(iter),liftz(iter))
if(rank == 0) then
- write(10,'(e14.5,e14.5)') current_time, drag(iter)
+! print *, 'drag: ', drag(iter)
+ write(10,'(e14.5,e14.5,e14.5,e14.5)') current_time, drag(iter), lifty(iter), liftz(iter)
end if
! Compute the new "penalized" vorticity
call curlDF4(velocity, vorticity,topo, grid_step)
! Computation of rhs of vorticity equation (on the grid) ...
call computeRHS(velocity, vorticity, rhs, topo, grid_step, nu)
- call MPI_BARRIER(MPI_COMM_WORLD,info)
+ !call MPI_BARRIER(MPI_COMM_WORLD,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,rhs, 3, info)
@@ -223,7 +237,7 @@ contains
! Remesh
call remesh(topo%ID, topo%mesh(1)%ID,domain_ghostsize,vorticity)
- if(mod(iter, 10).eq.0) then
+ if(mod(iter,100).eq.0) then
call printToVTK("post",iter, current_time,velocity,vorticity,topo,grid_step)
end if
current_time = current_time + time_step
@@ -234,7 +248,7 @@ contains
enddo
- close(10)
+ if(rank==0) close(10)
call ppm_finalize(info)
print *, 'end ppm simulation'
deallocate(drag)
diff --git a/HySoP/src/main/Particles.f90 b/HySoP/src/main/Particles.f90
index b9f9377f5..3c33f21bf 100755
--- a/HySoP/src/main/Particles.f90
+++ b/HySoP/src/main/Particles.f90
@@ -3,7 +3,7 @@ module Particles
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,rank
use mpi
use ppm_module_data, only : ppm_param_rmsh_kernel_mp4
use ppm_module_map_part
@@ -47,7 +47,7 @@ contains
! topo and mesh ids
integer, intent(in) :: topoid
integer, intent(in) :: meshid
- integer :: info, newNpart,rank,i
+ integer :: info, newNpart,i
! velocity on grid
real(mk), dimension(:,:,:,:,:), pointer :: vel
@@ -62,7 +62,6 @@ contains
npart = 0
last_npart = npart
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?
! --------> And deallocation will be handled by ppm also?
@@ -133,7 +132,7 @@ contains
! topo and mesh ids
integer, intent(in) :: topoid
integer, intent(in) :: meshid
- integer :: info, newNpart,rank
+ integer :: info, newNpart
integer, dimension(:), pointer :: ghostsize
! velocity on grid
real(mk), dimension(:,:,:,:,:), pointer :: vel
@@ -142,7 +141,6 @@ contains
info = 0
- call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info)
! Call ppm func to create the particles distribution
last_npart = npart
call ppm_rmsh_create_part(topoid,meshid,xp,npart,omp,dime,field_on_grid,cutoff,info,resetpos,&
@@ -152,24 +150,7 @@ contains
deallocate(rhsp,xp0)
allocate(rhsp(dime,npart),xp0(dime,npart))
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
!> --------> mesh to particles for velocity => done in ppm_rmsh_create_part ???
@@ -189,7 +170,7 @@ contains
! topo and mesh ids
integer, intent(in) :: topoid
integer, intent(in) :: meshid
- integer :: info!, rank
+ integer :: info
integer, dimension(:), pointer :: ghostsize
! velocity on grid
real(mk), dimension(:,:,:,:,:), pointer :: vel
@@ -197,8 +178,8 @@ contains
real(mk), dimension(:,:,:,:,:), pointer :: rhs
!> local loop indices
- integer :: i, newNpart,rank
- call MPI_COMM_RANK(MPI_COMM_WORLD, rank, info)
+ integer :: i, newNpart
+
do i=1,npart
xp0(1,i) = xp(1,i)
diff --git a/HySoP/src/main/Penalize.f90 b/HySoP/src/main/Penalize.f90
index 7543abd80..a1ebdd4da 100755
--- a/HySoP/src/main/Penalize.f90
+++ b/HySoP/src/main/Penalize.f90
@@ -7,7 +7,7 @@ module penalisation
private
- public :: penalisation_init, penalise_velocity,get_drag
+ public :: penalisation_init, penalise_velocity,get_dragAndLift
! Penalisation kernel function
real(mk), dimension(:,:,:,:), pointer :: chi => NULL()
@@ -20,7 +20,7 @@ module penalisation
real(mk) :: lambda
! Drag
- real(mk) :: local_drag
+ real(mk) :: local_drag,local_lifty,local_liftz
! Elementary volume
real(mk) :: dvol
@@ -105,6 +105,8 @@ contains
nullify(isublist, ndata)
local_drag = 0.0
+ local_lifty = 0.0
+ local_liftz = 0.0
dvol = meshStep(1)*meshStep(2)*meshStep(3)
end subroutine penalisation_init
@@ -123,7 +125,9 @@ contains
integer :: i,j,k, isub, nsublist,isubl
real(mk) :: coef
- local_drag = 0
+ local_drag = 0.0
+ local_lifty = 0.0
+ local_liftz = 0.0
nsublist = topo%nsublist
ndata => topo%mesh(1)%nnodes
@@ -137,6 +141,8 @@ contains
vel(2,i,j,k,isub)=vel(2,i,j,k,isub)*coef
vel(3,i,j,k,isub)=vel(3,i,j,k,isub)*coef
local_drag = local_drag + chi_sphere(i,j,k,isub)*vel(1,i,j,k,isub)*dvol
+ local_lifty = local_lifty + chi_sphere(i,j,k,isub)*vel(2,i,j,k,isub)*dvol
+ local_liftz = local_liftz + chi_sphere(i,j,k,isub)*vel(3,i,j,k,isub)*dvol
end do
end do
end do
@@ -146,17 +152,23 @@ contains
end subroutine penalise_velocity
-subroutine get_drag(drag)
+subroutine get_dragAndLift(drag,lifty,liftz)
real(mk), intent(out) :: drag
+ real(mk), intent(out) :: lifty
+ real(mk), intent(out) :: liftz
+
integer :: info
real(mk) :: coef,uinf
-! call MPI_BARRIER(MPI_COMM_WORLD,info)
uinf = 1.0
coef = 2./(uinf**2*2.*sphere_radius)
local_drag = coef*local_drag
- call MPI_AllReduce(local_drag, drag, 1, MPI_DOUBLE_PRECISION, MPI_SUM,MPI_COMM_WORLD,info)
+ local_lifty = coef*local_lifty
+ local_liftz = coef*local_liftz
+ call MPI_Reduce(local_drag, drag, 1, MPI_DOUBLE_PRECISION, MPI_SUM,0,MPI_COMM_WORLD,info)
+ call MPI_Reduce(local_lifty, lifty, 1, MPI_DOUBLE_PRECISION, MPI_SUM,0,MPI_COMM_WORLD,info)
+ call MPI_Reduce(local_liftz, liftz, 1, MPI_DOUBLE_PRECISION, MPI_SUM,0,MPI_COMM_WORLD,info)
-end subroutine get_drag
+end subroutine get_dragAndLift
end module Penalisation
diff --git a/HySoP/src/main/TestFunctions.f90 b/HySoP/src/main/TestFunctions.f90
index 224dbff48..30b5ffa35 100755
--- a/HySoP/src/main/TestFunctions.f90
+++ b/HySoP/src/main/TestFunctions.f90
@@ -1,6 +1,6 @@
module testsFunctions
- use client_data, only: mk
+ use client_data, only: mk,pi
use client_topology, only: ppm_t_topo
implicit none
@@ -137,9 +137,6 @@ contains
real(mk) :: x,y,z,cx,cy,cz,c2x,c2y,c2z,sx,sy,sz,s2x,s2y,s2z,x0,y0,z0
integer :: i,j,k,isub,isubl
- real(mk) :: pi
-
- pi = 4.0*atan(1.0_mk)
do isub=1,topo%nsublist
isubl=topo%isublist(isub)
diff --git a/HySoP/src/main/client_data.f90 b/HySoP/src/main/client_data.f90
index d641d930b..641424034 100755
--- a/HySoP/src/main/client_data.f90
+++ b/HySoP/src/main/client_data.f90
@@ -5,84 +5,6 @@ module client_data
integer, parameter :: mk = kind(1.0d0) ! double precision
integer, parameter :: dime = 3 ! Problem dimension
-
-!!$ !----------------------------------------------------------------------------!
-!!$ ! data particles
-!!$ !----------------------------------------------------------------------------!
-!!$ REAL(mk), DIMENSION(:,: ) , POINTER :: xp ! positions
-!!$ REAL(mk), DIMENSION(:,: ) , POINTER :: xp0 ! some odda positions
-!!$ REAL(mk), DIMENSION(:,: ) , POINTER :: wp ! vorticity
-!!$ REAL(mk), DIMENSION(:,: ) , POINTER :: wp0 ! some odda vorticity
-!!$ REAL(mk), DIMENSION(:,: ) , POINTER :: dwp ! vorticity
-!!$ REAL(mk), DIMENSION(:,: ) , POINTER :: up ! velocities
-
-
-!!$ !----------------------------------------------------------------------------!
-!!$ ! parallel stuff
-!!$ !----------------------------------------------------------------------------!
-!!$ CHARACTER(len=256) :: procname ! processor name
-!!$ INTEGER :: iprocname ! length of processor name
-!!$ INTEGER :: rank, nproc ! rank of this machine
-!!$ INTEGER :: ilogfile ! log file unit
-!!$ INTEGER :: debug ! debug flag for ppm
-!!$ INTEGER :: maxsubs, nsublist
-!!$
-!!$ CHARACTER(len=256) :: runtag
-!!$ INTEGER :: iruntag, comm, mpi_prec
-!!$ INTEGER :: coord(3),ndims(3)
-!!$
-!!$ !----------------------------------------------------------------------------!
-!!$ ! other things
-!!$ !----------------------------------------------------------------------------!
-!!$ REAL(mk), PARAMETER :: M_PI = 3.14159265358979
-!!$ REAL(mk) :: time, cutoff, tend
-!!$ INTEGER :: itime, krnl, ndump, itend
-!!$
-!!$ !----------------------------------------------------------------------------!
-!!$ ! time stepping
-!!$ !----------------------------------------------------------------------------!
-!!$ REAL(mk) :: dt
-!!$ REAL(mk) :: dt_max
-!!$ LOGICAL :: dt_adapt
-!!$ REAL(mk) :: lagcfl_omega, lagcfl_nablau, lagcfl_deform, fourier_diff
-!!$
-!!$ !----------------------------------------------------------------------------!
-!!$ ! physics
-!!$ !----------------------------------------------------------------------------!
-!!$ REAL(mk) :: nu
-!!$ REAL(mk) :: max_vorticity
-!!$ INTEGER :: g_istage
-!!$
-!!$ !----------------------------------------------------------------------------!
-!!$ ! parameters
-!!$ !----------------------------------------------------------------------------!
-!!$ INTEGER, PARAMETER :: vortex_prm_vorticity = 1
-!!$ INTEGER, PARAMETER :: vortex_prm_stream = 2
-!!$ INTEGER, PARAMETER :: vortex_prm_velocity = 3
-!!$ INTEGER, PARAMETER :: vortex_prm_dgammadt = 4
-!!$
-!!$ !----------------------------------------------------------------------------!
-!!$ ! Multigrid poisson solver
-!!$ !----------------------------------------------------------------------------!
-!!$ INTEGER :: maxlev
-!!$ !----------------------------------------------------------------------------!
-!!$ ! some other parameters
-!!$ !----------------------------------------------------------------------------!
-!!$ LOGICAL :: verbose
-!!$ REAL(mk) :: vortex_noise_amp
-!!$ REAL(mk) :: vortex_noise_basemode
-!!$ INTEGER :: vortex_noise_nmodes
-!!$
-!!$ !-----------------------------------------------------
-!!$ ! Trailing vortices
-!!$ !-----------------------------------------------------
-!!$ REAL(mk) :: trailvortex_b1, trailvortex_b2
-!!$ REAL(mk) :: trailvortex_a1, trailvortex_a2
-!!$ REAL(mk) :: trailvortex_z12
-!!$ REAL(mk) :: trailvortex_r, trailvortex_gamma
-!!$ LOGICAL :: trailvortex_noise_symmetry
-!!$ REAL(mk) :: trailvortex_noise_amp1, trailvortex_noise_amp2
-!!$ REAL(mk) :: trailvortex_noise_theta1, trailvortex_noise_theta2
-
-
+ real(mk), parameter :: pi = 4.0*atan(1.0_mk)
+ integer :: rank
END MODULE client_data
--
GitLab