forked from xcompact3d/Incompact3d
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathxcompact3d.f90
352 lines (274 loc) · 9.84 KB
/
xcompact3d.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
!Copyright (c) 2012-2022, Xcompact3d
!This file is part of Xcompact3d (xcompact3d.com)
!SPDX-License-Identifier: BSD 3-Clause
program xcompact3d
use var
use case
use transeq, only : calculate_transeq_rhs
use time_integrators, only : int_time
use navier, only : velocity_to_momentum, momentum_to_velocity, pre_correc, &
calc_divu_constraint, solve_poisson, cor_vel
use tools, only : restart, simu_stats, apply_spatial_filter, read_inflow
use turbine, only : compute_turbines
use ibm_param
use ibm, only : body
use genepsi, only : genepsi3d
use mhd, only : Bm,mhd_active,mhd_equation,test_magnetic, &
solve_poisson_mhd,mhd_sta
implicit none
call init_xcompact3d()
do itime=ifirst,ilast
!t=itime*dt
t=t0 + (itime0 + itime + 1 - ifirst)*dt
call simu_stats(2)
if (iturbine.ne.0) call compute_turbines(ux1, uy1, uz1)
if (iin.eq.3.and.mod(itime,ntimesteps)==1) then
call read_inflow(ux_inflow,uy_inflow,uz_inflow,itime/ntimesteps)
endif
if ((itype.eq.itype_abl.or.iturbine.ne.0).and.(ifilter.ne.0).and.(ilesmod.ne.0)) then
call filter(C_filter)
call apply_spatial_filter(ux1,uy1,uz1,phi1)
endif
do itr=1,iadvance_time
call set_fluid_properties(rho1,mu1)
call boundary_conditions(rho1,ux1,uy1,uz1,phi1,ep1)
if (imove.eq.1) then ! update epsi for moving objects
if ((iibm.eq.2).or.(iibm.eq.3)) then
call genepsi3d(ep1)
else if (iibm.eq.1) then
call body(ux1,uy1,uz1,ep1)
endif
endif
call calculate_transeq_rhs(drho1,dux1,duy1,duz1,dphi1,rho1,ux1,uy1,uz1,ep1,phi1,divu3)
#ifdef DEBG
call check_transients()
#endif
if (ilmn) then
!! XXX N.B. from this point, X-pencil velocity arrays contain momentum (LMN only).
call velocity_to_momentum(rho1,ux1,uy1,uz1)
endif
call int_time(rho1,ux1,uy1,uz1,phi1,drho1,dux1,duy1,duz1,dphi1)
call pre_correc(ux1,uy1,uz1,ep1)
call calc_divu_constraint(divu3,rho1,phi1)
call solve_poisson(pp3,px1,py1,pz1,rho1,ux1,uy1,uz1,ep1,drho1,divu3)
call cor_vel(ux1,uy1,uz1,px1,py1,pz1)
if(mhd_active .and. mhd_equation) then
call solve_poisson_mhd()
endif
if (ilmn) then
call momentum_to_velocity(rho1,ux1,uy1,uz1)
!! XXX N.B. from this point, X-pencil velocity arrays contain velocity (LMN only).
!! Note - all other solvers work on velocity always
endif
call test_flow(rho1,ux1,uy1,uz1,phi1,ep1,drho1,divu3)
if(mhd_active) call test_magnetic
enddo !! End sub timesteps
if(mhd_active) call mhd_sta(ux1,uy1,uz1)
!
call restart(ux1,uy1,uz1,dux1,duy1,duz1,ep1,pp3(:,:,:,1),phi1,dphi1,px1,py1,pz1,rho1,drho1,mu1,1)
call simu_stats(3)
call postprocessing(rho1,ux1,uy1,uz1,pp3,phi1,ep1)
enddo !! End time loop
call finalise_xcompact3d()
end program xcompact3d
!########################################################################
!########################################################################
subroutine init_xcompact3d()
use MPI
use decomp_2d
use decomp_2d_io, only : decomp_2d_io_init
USE decomp_2d_poisson, ONLY : decomp_2d_poisson_init
use case
use sandbox, only : init_sandbox
use forces
use var
use navier, only : calc_divu_constraint
use tools, only : test_speed_min_max, test_scalar_min_max, &
restart, &
simu_stats, compute_cfldiff, &
init_inflow_outflow, read_inflow
use param, only : ilesmod, jles,itype
use param, only : irestart
use variables, only : nx, ny, nz, nxm, nym, nzm
use variables, only : p_row, p_col
use variables, only : nstat, nvisu, ilist
use les, only: init_explicit_les
use turbine, only: init_turbines
use visu, only : visu_init, visu_ready
use genepsi, only : genepsi3d, epsi_init
use ibm, only : body
use probes, only : init_probes
use mhd, only: mhd_active,mhd_init
implicit none
integer :: ierr
integer :: nargin, FNLength, status, DecInd
logical :: back
character(len=80) :: InputFN, FNBase
!! Initialise MPI
call MPI_INIT(ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD,nrank,ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,ierr)
! Handle input file like a boss -- GD
nargin=command_argument_count()
if (nargin <1) then
InputFN='input.i3d'
if (nrank==0) write(*,*) 'Xcompact3d is run with the default file -->', trim(InputFN)
elseif (nargin >= 1) then
call get_command_argument(1,InputFN,FNLength,status)
back=.true.
FNBase=inputFN((index(InputFN,'/',back)+1):len(InputFN))
DecInd=index(FNBase,'.',back)
if (DecInd >1) then
FNBase=FNBase(1:(DecInd-1))
end if
if (nrank==0) write(*,*) 'Xcompact3d is run with the provided file -->', trim(InputFN)
endif
#ifdef ADIOS2
if (nrank .eq. 0) then
print *, " WARNING === WARNING === WARNING === WARNING === WARNING"
print *, " WARNING: Running Xcompact3d with ADIOS2"
print *, " this is currently experimental"
print *, " for safety of results it is recommended"
print *, " to run the default build as this feature"
print *, " is developed. Thank you for trying it."
print *, " WARNING === WARNING === WARNING === WARNING === WARNING"
endif
#endif
call parameter(InputFN)
call decomp_2d_init(nx,ny,nz,p_row,p_col)
call decomp_2d_io_init()
call init_coarser_mesh_statS(nstat,nstat,nstat,.true.) !start from 1 == true
call init_coarser_mesh_statV(nvisu,nvisu,nvisu,.true.) !start from 1 == true
!div: nx ny nz --> nxm ny nz --> nxm nym nz --> nxm nym nzm
call decomp_info_init(nxm, nym, nzm, ph1)
call decomp_info_init(nxm, ny, nz, ph4)
!gradp: nxm nym nzm -> nxm nym nz --> nxm ny nz --> nx ny nz
call decomp_info_init(nxm, ny, nz, ph2)
call decomp_info_init(nxm, nym, nz, ph3)
call init_variables()
call schemes()
call decomp_2d_poisson_init()
call decomp_info_init(nxm,nym,nzm,phG)
if (ilesmod.ne.0) then
if (jles.gt.0) call init_explicit_les()
endif
if ((iibm.eq.2).or.(iibm.eq.3)) then
call genepsi3d(ep1)
else if (iibm.eq.1) then
call epsi_init(ep1)
call body(ux1,uy1,uz1,ep1)
endif
if (iforces.eq.1) then
call init_forces()
if (irestart==1) then
call restart_forces(0)
endif
endif
!####################################################################
! initialise mhd
if (mhd_active) call mhd_init()
!####################################################################
! initialise visu
if (ivisu.ne.0) then
call visu_init()
call visu_case_init() !! XXX: If you get error about uninitialised IO, look here.
!! Ensures additional case-specific variables declared for IO
call visu_ready()
end if
! compute diffusion number of simulation
call compute_cfldiff()
!####################################################################
if (irestart==0) then
call init(rho1,ux1,uy1,uz1,ep1,phi1,drho1,dux1,duy1,duz1,dphi1,pp3,px1,py1,pz1)
itime = 0
call preprocessing(rho1,ux1,uy1,uz1,pp3,phi1,ep1)
else
itr=1
if (itype == itype_sandbox) then
call init_sandbox(ux1,uy1,uz1,ep1,phi1,1)
end if
call restart(ux1,uy1,uz1,dux1,duy1,duz1,ep1,pp3(:,:,:,1),phi1,dphi1,px1,py1,pz1,rho1,drho1,mu1,0)
endif
if ((ioutflow.eq.1).or.(iin.eq.3)) then
call init_inflow_outflow()
if ((irestart==1).and.(iin.eq.3)) then
call read_inflow(ux_inflow,uy_inflow,uz_inflow,itime/ntimesteps)
endif
end if
if ((iibm.eq.2).or.(iibm.eq.3)) then
call genepsi3d(ep1)
else if (iibm.eq.1) then
call body(ux1,uy1,uz1,ep1)
endif
if (mod(itime, ilist) == 0 .or. itime == ifirst) then
call test_speed_min_max(ux1,uy1,uz1)
if (iscalar==1) call test_scalar_min_max(phi1)
endif
call simu_stats(1)
call calc_divu_constraint(divu3, rho1, phi1)
call init_probes()
if (iturbine.ne.0) call init_turbines(ux1, uy1, uz1)
if (itype==2) then
if(nrank.eq.0)then
open(42,file='time_evol.dat',form='formatted')
endif
endif
if (itype==5) then
if(nrank.eq.0)then
open(38,file='forces.dat',form='formatted')
endif
endif
endsubroutine init_xcompact3d
!########################################################################
!########################################################################
subroutine finalise_xcompact3d()
use MPI
use decomp_2d_mpi
use decomp_2d
use decomp_2d_io, only : decomp_2d_io_finalise
use tools, only : simu_stats
use param, only : itype, jles, ilesmod
use probes, only : finalize_probes
use visu, only : visu_finalise
use les, only: finalise_explicit_les
use mhd, only: mhd_active, mhd_fin
implicit none
integer :: ierr
if (itype==2) then
if(nrank.eq.0)then
close(42)
endif
endif
if (itype==5) then
if(nrank.eq.0)then
close(38)
endif
endif
call simu_stats(4)
call finalize_probes()
call visu_finalise()
if (mhd_active) call mhd_fin()
if (ilesmod.ne.0) then
if (jles.gt.0) call finalise_explicit_les()
endif
call decomp_2d_io_finalise()
call decomp_2d_finalize
CALL MPI_FINALIZE(ierr)
endsubroutine finalise_xcompact3d
subroutine check_transients()
use decomp_2d_constants, only : mytype
use mpi
use var
implicit none
real(mytype) :: dep, dep1
integer :: code
dep=maxval(abs(dux1))
call MPI_ALLREDUCE(dep,dep1,1,real_type,MPI_MAX,MPI_COMM_WORLD,code)
if (nrank == 0) write(*,*)'## MAX dux1 ', dep1
dep=maxval(abs(duy1))
call MPI_ALLREDUCE(dep,dep1,1,real_type,MPI_MAX,MPI_COMM_WORLD,code)
if (nrank == 0) write(*,*)'## MAX duy1 ', dep1
dep=maxval(abs(duz1))
call MPI_ALLREDUCE(dep,dep1,1,real_type,MPI_MAX,MPI_COMM_WORLD,code)
if (nrank == 0) write(*,*)'## MAX duz1 ', dep1
end subroutine check_transients