From ebef18826f69aac8903c09905e58137185d53324 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Fri, 23 Sep 2011 07:58:43 +0000
Subject: [PATCH] Last (minor) updates in NS
---
HySoP/src/main/NavierStokes3D.f90 | 40 +++++++++++++------------------
HySoP/src/main/Penalize.f90 | 6 +++--
2 files changed, 21 insertions(+), 25 deletions(-)
diff --git a/HySoP/src/main/NavierStokes3D.f90 b/HySoP/src/main/NavierStokes3D.f90
index d4d2410f7..46f1c9412 100755
--- a/HySoP/src/main/NavierStokes3D.f90
+++ b/HySoP/src/main/NavierStokes3D.f90
@@ -126,9 +126,9 @@ contains
real(mk) ::Re
real(mk) ::mem_used,t1,t2
- Re = 210.
+ Re = 133.4
initial_time = 0.0
- final_time = 70.
+ final_time = 60.
! Initialisation stuff
t1 = MPI_WTIME()
@@ -151,11 +151,12 @@ contains
! allocate(err1(dime,topo%nsublist),err2(dime,topo%nsublist),errinf(dime,topo%nsublist))
print *, '[',rank,'] space step: ', grid_step
print *, '[',rank,'] mesh dim: ', topo%mesh(1)%nnodes
+ if(rank==0) print *, 'Reynolds : ', Re
mem_used=3*sizeof(velocity)*1d-6
print * , rank, ']] Theoretical used memory (MB) for fields : ', mem_used
- ! do while(final_time>initial_time)
+ ! do while(final_time>initial_time)
!enddo
! Time loop
call timeLoop(initial_time,final_time,Re)
@@ -177,7 +178,7 @@ contains
real(mk) :: current_time,elapsed_time,time_step
integer :: iter,max_iter
logical :: resetpos
- real(mk) :: fit_velo,nu, timer
+ real(mk) :: fit_velo,nu
! real(mk), dimension(:,:), pointer :: err1, err2, errinf
real(mk), dimension(dime) :: diagnostics ! (1): drag, (2):lifty, (3):liftz
@@ -186,15 +187,15 @@ contains
iter = 1
! offset for velocity to set a desired average value for flow rate.
fit_velo = (boundary_layer_max(3)-boundary_layer_min(3))/lengths(3)
- ! Allocate drag and lifts vectors (one value per time step)
- ! max_iter= int((final_time -initial_time)/time_step) + 1
- ! allocate(drag(max_iter),lifty(max_iter),liftz(max_iter))
+
if(rank==0) open(10,file='diagnostics') ! Output only for one process
! Synchronise all procs after initialization process ...
call MPI_BARRIER(MPI_COMM_WORLD,info)
print *, '[',rank,'] -------------------------------> start simulation'
+ ! Initial time step
time_step=0.005
+
resetpos = .TRUE.
do while(current_time <= final_time)
@@ -208,15 +209,13 @@ contains
! Solve poisson to find velocity
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 MPI_BARRIER(MPI_COMM_WORLD,info)
-
velocity(1,:,:,:,:) = velocity(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
+!!$ if( current_time>3.0 .and. current_time<4.0) then
+!!$ call perturb(velocity, current_time)
+!!$ end if
! Penalize velocity and map result to all procs
call penalise_velocity(velocity,topo,time_step)
@@ -249,16 +248,9 @@ contains
! Initialize the particles according to the last computed vorticity, velocity and rhs
call update_particles(vorticity,resetpos,topo%ID,topo%mesh(1)%ID,velocity)
- ! Try to dealloc/alloc vorticity to lower memory use.
- ! Failed during remesh call. It seems that ppm needs the same adress for vorticity in update and remesh.
- ! deallocate(vorticity)
-
! Integrate
call push_particles(time_step,topo%ID,topo%mesh(1)%ID,domain_ghostsize,velocity,rhs)
- timer = MPI_Wtime()
- ! allocate(vorticity(dime,size(velocity,2),size(velocity,3),size(velocity,4),size(velocity,5)))
-
! Remesh
call remesh(topo%ID, topo%mesh(1)%ID,domain_ghostsize,vorticity)
@@ -267,16 +259,18 @@ contains
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)
-
+
+ ! Update time
current_time = current_time + time_step
! Compute the new time step according to vorticity maximum value
call updateTimeStep(time_step,vorticity)
! Output every 10 time unit ...
- if(mod(current_time,10.).lt.time_step) then
- call printToVTK("run",iter,current_time,velocity,vorticity,topo,grid_step)
- end if
+ !if(mod(current_time,10.).lt.time_step) then
+ ! call printToVTK("run",iter,current_time,velocity,vorticity,topo,grid_step)
+ !end if
+
!current_time = final_time
iter = iter+1
print *, "[", rank,"] time : ", MPI_WTime()-elapsed_time
diff --git a/HySoP/src/main/Penalize.f90 b/HySoP/src/main/Penalize.f90
index 2fe26c911..4f99e073e 100755
--- a/HySoP/src/main/Penalize.f90
+++ b/HySoP/src/main/Penalize.f90
@@ -60,8 +60,10 @@ contains
real(mk) :: dist
integer :: sizeMaxChi,count_sphere,count_boundary
! Set lambda value
- lambda = 6.d5
-
+
+ lambda = 6.e6
+
+ if(rank==0) print *, 'Lambda for penalization: ', lambda
! Set sphere position and radius
sphere_radius = 0.5
sphere_pos = (boundary_layer_min + boundary_layer_max)*0.5
--
GitLab