-
Notifications
You must be signed in to change notification settings - Fork 12
/
main.f90
356 lines (283 loc) · 11.3 KB
/
main.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
program x_code
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
use m_stats
use m_timing
use m_force
use m_rand_knuth
use m_particles
use m_filter_xfftw
use m_les
implicit none
integer :: n
character :: sym
call m_timing_init ! Setting the time zero
call m_openmpi_init
call m_io_init
call m_parameters_init
call m_les_init
call m_fields_init
call m_work_init
! allocating and initializing FFTW arrays
call x_fftw_allocate(1)
call x_fftw_init
call m_stats_init
call m_force_init
! allocating and initializing particles
if (task.eq.'parts') then
call particles_init
end if
write(out,*) "IN THE PROGRAM."
call flush(out)
! initializing the random number generator
! call rand_knuth_init
! getting the wallclock runlimit for the job
call get_job_runlimit
!-----------------------------------------------------------------------
! Starting from the beginning or from the saved flowfield
!-----------------------------------------------------------------------
if(ITMIN.eq.0) then
call begin_new
else
call begin_restart
endif
! Initializing the LES stuff
if (les) call m_les_begin
! checking divergence
if (task.eq.'hydro') call divergence
! indicators whether to use first-order in time
! for velocities and scalars
fov = .true.
fos = .true.
! need to dealias the fields at the beginning
if (task.eq.'hydro') call dealias_all
!********************************************************************************
! call benchmark
!********************************************************************************
!================================================================================
! MAIN CYCLE
!================================================================================
do 100 ITIME=ITMIN+1,ITMAX
! getting the file extension for current iteration
call get_file_ext
!--------------------------------------------------------------------------------
! now performing the core of the cycle.
! This is done with "if" rather than "select case" because if we're not
! splitting tasks then we want everything to be done consequently by the
! same set of processors.
!
! All the syncronization calls (fields_to_parts, fields_to_stats) will be
! called only if (task_split).
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
! HYDRO PART
! note that even in the case where there is no task splitting,
! 'hydro' part is still there. all processors will have task = 'hydro'
!--------------------------------------------------------------------------------
hydro: if (task.eq.'hydro') then
! ------------------------------------------------------------
! taking care of rescaling when running decaying turbulence
! if the time just was divisible by TRESCALE
! ------------------------------------------------------------
if (flow_type.eq.0 .and. floor((time-dt)/TRESCALE) .lt. floor(time/TRESCALE)) then
! ...and if we haven't rescaled NRESCALE times
if (floor(time/TRESCALE) .le. NRESCALE .and. itime.ne.1) then
write(out,*) "MAIN: Rescaling velocities"
call flush(out)
call velocity_rescale
! after rescaling, the time-sceping needs to be first order
fov = .true.; fos = .true.
if (.not. task_split .and. mod(itime,iprint1).eq.0) call stat_main
end if
end if
! RHS for scalars
call rhs_scalars
! now the velocities in x-space are contained in wrk1...3
! if we are moving particles, then we want to send the velocity field
! to the "parts" part of the code
if (task_split) call fields_to_parts
! advance scalars - either Euler or Adams-Bashforth
if (int_scalars .or. n_les > 0) then
call flush(out)
n = 3 + n_scalars + n_les
if (fos) then
rhs_old(:,:,:,4:n) = wrk(:,:,:,4:n)
fields(:,:,:,4:n) = fields(:,:,:,4:n) + dt * rhs_old(:,:,:,4:n)
fos = .false.
else
fields(:,:,:,4:n) = fields(:,:,:,4:n) + &
dt * ( 1.5d0 * wrk(:,:,:,4:n) - 0.5d0 * rhs_old(:,:,:,4:n) )
rhs_old(:,:,:,4:n) = wrk(:,:,:,4:n)
end if
end if
! RHS for velocities
call rhs_velocity
! adding forcing, if computing a forced flow
if (flow_type.eq.1) call force_velocity
! advance velocity - either Euler or Adams-Bashforth
if (fov) then
rhs_old(:,:,:,1:3) = wrk(:,:,:,1:3)
fields(:,:,:,1:3) = fields(:,:,:,1:3) + dt * rhs_old(:,:,:,1:3)
fov = .false.
else
fields(:,:,:,1:3) = fields(:,:,:,1:3) + &
dt * ( 1.5d0 * wrk(:,:,:,1:3) - 0.5d0 * rhs_old(:,:,:,1:3) )
rhs_old(:,:,:,1:3) = wrk(:,:,:,1:3)
end if
! solve for pressure and update velocities so they are incompressible
call pressure
! advance the time
TIME = TIME + DT
! write the restart file if it's the time
if (mod(itime,IPRINT2).eq.0) call restart_write_parallel
! change the timestep in case we're running with variable timestep
if (variable_dt) call my_dt
! CPU usage statistics
if (mod(itime,iprint1).eq.0) then
call m_timing_check
if (mod(itime,iwrite4).eq.0) then
sym = "*"
else
sym = " "
end if
write(out,9000) itime,time,dt,courant,cpu_hrs,cpu_min,cpu_sec,&
sym,les_model_name
call flush(out)
end if
if (mod(itime,iprint1).eq.0 .or. mod(itime,iwrite4).eq.0) then
! send the velocities to the "stats" part of the code for statistics
if (task_split) call fields_to_stats
! checking if we need to stop the calculations due to simulation time
if (TIME.gt.TMAX) call my_exit(1)
! checking if we need to start advancing scalars
if (n_scalars.gt.0 .and. .not.int_scalars .and. time.gt.TSCALAR) then
int_scalars = .true.
call init_scalars
write(out,*) "Starting to move the scalars."
call flush(out)
end if
end if
end if hydro
!--------------------------------------------------------------------------------
! STATISTICS PART
!--------------------------------------------------------------------------------
stats: if (task.eq.'stats' .or. .not.task_split) then
if (mod(itime,iprint1).eq.0 .or. mod(itime,iwrite4).eq.0) then
! if this is a separate set of processors, then...
stats_task_split: if (task_split) then
! checking if we need to stop the calculations due to simulation time
if (TIME.gt.TMAX) call my_exit(1)
end if stats_task_split
! these are executed regardless of the processor configuration
if (task_split) call fields_to_stats
if (mod(itime,iprint1).eq.0) call stat_main
if (mod(itime,iwrite4).eq.0) call io_write_4
end if
end if stats
!--------------------------------------------------------------------------------
! PARTICLE PARTS
! NOTE: This is not enabled to work when not task_split.
! Need to return to it later.
! Currently the particles can be calculated only if we split the tasks due to
! requirements on the wrk array sizes in the particle interpolation routines.
!--------------------------------------------------------------------------------
particles: if (task.eq.'parts') then
call fields_to_parts
if (int_particles) then
call particles_move
if (mod(itime,iwrite4).eq.0) call particles_restart_write_binary
end if
if (mod(itime,iprint1).eq.0 .or. mod(itime,iwrite4).eq.0) then
if (TIME.gt.TMAX) call my_exit(1)
end if
end if particles
!!$!--------------------------------------------------------------------------------
!!$! OTHER PARTS
!!$!--------------------------------------------------------------------------------
!!$
!!$
!!$ write(out,*) "skipping the time step",ITIME
!!$ call flush(out)
!!$
!!$
!--------------------------------------------------------------------------------
! COMMON PARTS
!--------------------------------------------------------------------------------
! every 10 iterations checking
! 1) for the run time: are we getting close to the job_runlimit?
! 2) for the user termination: is there a file "stop" in directory?
if (mod(ITIME,10).eq.0) then
! synchronize all processors, hard
!!$ call MPI_BARRIER(MPI_COMM_WORLD,mpi_err)
if (myid_world.eq.0) call m_timing_check
count = 1
call MPI_BCAST(cpu_min_total,count,MPI_INTEGER4,0,MPI_COMM_WORLD,mpi_err)
! allowing 5 extra minutes for writing restart file
! note that for large-scale calculations (e.g. 1024^3)
! the restart writing time can be long (up to 20 minutes or so).
! this should be taken care of in the job submission script
! via file job_parameters.txt
if (cpu_min_total+5 .gt. job_runlimit) call my_exit(2)
! user termination. If the file "stop" is in the directory, stop
inquire(file='stop',exist=there)
if (there) call my_exit(3)
end if
100 continue
!================================================================================
!--------------------------------------------------------------------------------
! In a case when we've gone to ITMAX, write the restart file
!--------------------------------------------------------------------------------
ITIME = ITIME-1
if (task.eq.'hydro') call restart_write_parallel
call my_exit(0)
call m_openmpi_exit
stop
9000 format('ITIME=',i6,3x,'TIME=',f8.4,4x,'DT=',f8.5,3x,'Courn= ',f6.4, &
2x,'CPU:(',i4.4,':',i2.2,':',i2.2,')',x,a1,x,a3)
end program x_code
!=============================================================================
subroutine benchmark
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
use m_stats
use m_timing
use m_force
use m_rand_knuth
use m_particles
use m_filter_xfftw
use m_les
implicit none
integer :: n, nmax
call m_timing_init
benchmarking = .true.
bm = 0
wrk(:,:,:,1) = fields(:,:,:,1)
do n = 1, nmax
call xfft3d(1,1)
call xfft3d(-1,1)
end do
if (myid.eq.0) then
write(out,*) "BENF: statistics on forward transform"
write(out,*) "BENF: R2C: ", bm(1)/nmax
write(out,*) "BENF: T13: ", bm(2)/nmax
write(out,*) "BENF: C2C: ", bm(3)/nmax
write(out,*) "BENF: ====="
write(out,*) "BENF: TOT: ", bm(11)/nmax
write(out,*) "BENB: statistics on backward transform"
write(out,*) "BENB: C2C: ", bm(4)/nmax
write(out,*) "BENB: T13: ", bm(5)/nmax
write(out,*) "BENB: C2R: ", bm(6)/nmax
write(out,*) "BENB: ====="
write(out,*) "BENB: TOT: ", bm(12)/nmax
end if
close(out)
stop
end subroutine benchmark