pyCALC-LES
Loading...
Searching...
No Matches
exec-pyCALC-LES-GPU.py
Go to the documentation of this file.
1gpu = False
2import numpy
3from cupyx.scipy import sparse
4import sys
5import time
6import pyamg
7from cupyx.scipy.sparse import spdiags,linalg,eye
8
10 global acrank,acrank_conv, acrank_conv_keps, acrank_conv_kom,amg_cycle, amg_cycle_phi, amg_relax, amg_relax_phi, \
11 blend,cdes,c_eps,c_eps_1,c_eps_2, c_l,c_omega_1, c_omega_2, cmu,c_t, coeff_v, coeff_w,\
12 c_omega_1_sst_1, c_omega_2_sst_1, c_omega_1_sst_2,c_omega_2_sst_2, \
13 convergence_limit_eps, convergence_limit_k, convergence_limit_om, convergence_limit_p, convergence_limit_t, convergence_limit_u, \
14 convergence_limit_v, convergence_limit_w,cr_sst, convergence_limit_gpu, \
15 cyclic_x, cyclic_z, dt, dist3d,dz, dmin_synt,embedded,eps_bc_east,eps_bc_east_type,eps_bc_north,eps_bc_north_type, eps_bc_south, \
16 eps_bc_south_type, eps_bc_west, eps_bc_west_type, eps_bc_low, eps_bc_high, eps_bc_low_type,eps_bc_high_type,eps_min,fkmin_limit,\
17 fx, fy,fz,imon,itstep_stats, gpu, save_average_z, i_s_fair, k_s_fair, i_fair,i_fair_embedd, \
18 i_block_start,i_block_end,j_block_start,j_block_end,\
19 itstep_save,itstep_start,jl0,jmirror_synt,jmon,kappa,k_bc_east,k_bc_east_type,k_bc_north,k_bc_north_type,k_bc_south,\
20 k_bc_south_type,k_bc_west,k_bc_west_type,k_bc_low,k_bc_high,k_bc_low_type,k_bc_high_type,keps,keps_des,k_eq_les,kmon,kom, k_min,\
21 kom_peng,kom_des,launder,L_t_synt,\
22 maxit,min_iter, norm_order, ni,nj,nk,nmodes_synt,nsweep_keps,nsweep_kom, nsweep_t, \
23 nsweep_vel, ntstep, om_bc_east, om_bc_east_type, om_bc_north, om_bc_north_type, \
24 om_bc_south, om_bc_south_type, om_bc_west, om_bc_west_type, om_bc_low,om_bc_high, om_bc_low_type, om_bc_high_type, om_min,\
25 p_bc_east, p_bc_east_type, p_bc_north, p_bc_north_type, p_bc_south, p_bc_south_type, p_bc_west, p_bc_west_type, p_bc_low,\
26 p_bc_high, p_bc_low_type, p_bc_high_type,pans, prand_lam,prand_eps,\
27 prand_k_sst_1, prand_k_sst_2, prand_omega_sst_1, prand_omega_sst_2, prand_k,prand_omega,prand_t, \
28 resnorm_p,resnorm_vel,resnorm_t,restart,save,save_vtk_movie,scheme,scheme_turb,scheme_t, smag,solver_vel, solver_p, \
29 solver_t,t_bc_east, t_bc_east_type, t_bc_north, t_bc_north_type, t_bc_sotth, t_bc_south_type, t_bc_west, t_bc_west_type, \
30 t_bc_low,t_bc_high, t_bc_low_type, t_bc_high_type, temp, \
31 solver_turb,solverx,sormax, s2,sst,sst_sst_uv_limit,u_bc_east, u_bc_east_type, u_bc_north, u_bc_north_type, \
32 u_bc_south, u_bc_south_type, u_bc_west, u_bc_west_type, \
33 u_bc_low,u_bc_high,u_bc_low_type,u_bc_high_type,urfvis, v_bc_east, v_bc_east_type, v_bc_north, v_bc_north_type, v_bc_south, \
34 v_bc_south_type,v_bc_west,v_bc_west_type,v_bc_low,v_bc_high,v_bc_low_type,v_bc_high_type,viscos, vol,vtk,vtk_save,\
35 vtk_file_name,w_bc_east,w_bc_east_type,w_bc_north,w_bc_north_type,\
36 w_bc_south, w_bc_south_type, w_bc_west, w_bc_west_type, w_bc_low, w_bc_high, w_bc_low_type, w_bc_high_type,\
37 wale, x_embed, x2d, xp2d, y2d, yp2d, z,zp,zmax
38
39 import numpy as np
40 import sys
41
42# N.B. All variables that are set in this module must be included in the 'return' statement at the last line
43
44 gpu = False
45 scheme='h' #hybrid
46 scheme_turb='u' #hybrid upwind-central
47 acrank=1.0 # for pressure gradient
48 acrank_conv=1.0 # for convection-diffusion
49 acrank_conv_kom=1 # for convection-diffusion
50 acrank_conv_keps=1 # for convection-diffusion
51# scheme_turb='h' #hybrid upwind-central
52
53
54 cmu=0.09
55 pans = False
56 keps = False
57 kom_des = False
58 keps_des = False
59 kom = True
60 wale = False
61 smag = False
62 c_eps_1=1.5
63 c_eps_2=1.9
64 cmu=0.09
65 c_omega_1= 5./9.
66 c_omega_2=3./40.
67 prand_omega=2.0
68 prand_eps=1.4
69 prand_k=1.4
70 jl0=0
71
72 if keps:
73 prand_k=1.4
74
75 if pans:
76 prand_k=-1.4 # will be multiplied by fk3d in coeff()
77 prand_eps=-1.4 # will be multiplied by fk3d in coeff()
78
79 if kom or kom_des:
80 prand_k=2.0
81 if smag:
82 cmu=0.1
83
84 restart = False
85 save = True
86
87 viscos=1/5200
88
89 urfvis=0.5
90
91 maxit=5
92 min_iter=1
93 sormax=1e-3
94
95 amg_relax='default'
96 solver_vel='lgmres'
97 solver_turb='lgmres'
98 solver_p='pyamg'
99 amg_relax='fgmres'
100 amg_cycle='F'
101 nsweep_vel=50
102 nsweep_keps=50
103 nsweep_kom=50
104 convergence_limit_u=1e-5
105 convergence_limit_v=1e-5
106 convergence_limit_w=1e-5
107 convergence_limit_eps=1e-4
108 convergence_limit_k=-1e-4
109 convergence_limit_om=1e-8
110 convergence_limit_om=-1e-5
111 convergence_limit_p=5e-4
112
113 imon=0
114 jmon=0
115 kmon=0
116
117 ntstep=1000
118 uin=20
119 dt=0.5*(x2d[1,0]-x2d[0,0])*xp.ones(ntstep)/uin
120 dt=4*(x2d[1,0]-x2d[0,0])*xp.ones(ntstep)/uin
121 itstep_start=ntstep-100
122 itstep_save=2000 # save every itstep_save timestep
123 itstep_stats=1 # time average every itstep_stats timestep
124 vtk=False
125
126
127 resnorm_p=uin*zmax*y2d[1,-1]
128 resnorm_vel=uin**2*zmax*y2d[1,-1]
129
130
131 cyclic_x = True
132 cyclic_z = False
133 cyclic_z = True
134
135# synthetic inlet fluct
136 L_t_synt=0.2
137 nmodes_synt=150
138 jmirror_synt=int(nj/2) # mirror vsynt at node jmirror; jmirror=0 means no mirroring
139 dmin_synt=dz/2
140
141# boundary conditions for u
142 u_bc_west=xp.zeros((nj,nk))
143 u_bc_east=xp.zeros((nj,nk))
144 u_bc_south=xp.zeros((ni,nk))
145 u_bc_north=xp.zeros((ni,nk))
146 u_bc_low=xp.zeros((ni,nj))
147 u_bc_high=xp.zeros((ni,nj))
148
149 u_bc_west_type='d'
150 u_bc_east_type='n'
151 u_bc_south_type='d'
152 u_bc_north_type='d'
153 u_bc_low_type='d'
154 u_bc_high_type='d'
155
156
157# boundary conditions for v
158 v_bc_west=xp.zeros((nj,nk))
159 v_bc_east=xp.zeros((nj,nk))
160 v_bc_south=xp.zeros((ni,nk))
161 v_bc_north=xp.zeros((ni,nk))
162 v_bc_low=xp.zeros((ni,nj))
163 v_bc_high=xp.zeros((ni,nj))
164
165 v_bc_west_type='d'
166 v_bc_east_type='n'
167 v_bc_south_type='d'
168 v_bc_north_type='d'
169 v_bc_low_type='d'
170 v_bc_high_type='d'
171
172# boundary conditions for w
173 w_bc_west=xp.zeros((nj,nk))
174 w_bc_east=xp.zeros((nj,nk))
175 w_bc_south=xp.zeros((ni,nk))
176 w_bc_north=xp.zeros((ni,nk))
177 w_bc_low=xp.zeros((ni,nj))
178 w_bc_high=xp.zeros((ni,nj))
179
180 w_bc_west_type='d'
181 w_bc_east_type='n'
182 w_bc_south_type='d'
183 w_bc_north_type='d'
184 w_bc_low_type='d'
185 w_bc_high_type='d'
186
187# boundary conditions for p
188 p_bc_west=xp.zeros((nj,nk))
189 p_bc_east=xp.zeros((nj,nk))
190 p_bc_south=xp.zeros((ni,nk))
191 p_bc_north=xp.zeros((ni,nk))
192 p_bc_low=xp.zeros((ni,nj))
193 p_bc_high=xp.zeros((ni,nj))
194
195 p_bc_west_type='n'
196 p_bc_east_type='n'
197 p_bc_south_type='n'
198 p_bc_north_type='n'
199 p_bc_low_type='n'
200 p_bc_high_type='n'
201
202# boundary conditions for k
203 k_bc_west=xp.zeros((nj,nk))
204 k_bc_east=xp.zeros((nj,nk))
205 k_bc_south=xp.zeros((ni,nk))
206 k_bc_north=xp.zeros((ni,nk))
207 k_bc_low=xp.zeros((ni,nj))
208 k_bc_high=xp.zeros((ni,nj))
209
210 k_bc_west_type='d'
211 k_bc_east_type='d'
212 k_bc_south_type='d'
213 k_bc_north_type='d'
214 k_bc_low_type='d'
215 k_bc_high_type='d'
216
217# boundary conditions for eps
218 eps_bc_west=xp.zeros((nj,nk))
219 eps_bc_east=xp.zeros((nj,nk))
220 eps_bc_south=xp.zeros((ni,nk))
221 eps_bc_north=xp.zeros((ni,nk))
222 eps_bc_low=xp.zeros((ni,nj))
223 eps_bc_high=xp.zeros((ni,nj))
224
225 eps_bc_west_type='d'
226 eps_bc_east_type='d'
227 eps_bc_south_type='d'
228 eps_bc_north_type='d'
229 eps_bc_low_type='n'
230 eps_bc_high_type='n'
231
232# boundary conditions for omega
233 om_bc_west=xp.zeros((nj,nk))
234 om_bc_east=xp.zeros((nj,nk))
235 om_bc_south=xp.zeros((ni,nk))
236 om_bc_north=xp.zeros((ni,nk))
237 om_bc_low=xp.zeros((ni,nj))
238 om_bc_high=xp.zeros((ni,nj))
239
240 xwall_s=0.5*(x2d[0:-1,0]+x2d[1:,0])
241 ywall_s=0.5*(y2d[0:-1,0]+y2d[1:,0])
242 dist2_s=(yp2d[:,0]-ywall_s)**2+(xp2d[:,0]-xwall_s)**2
243 om_bc_south=6*viscos/c_omega_2/dist2_s
244
245# make it 2D
246 om_bc_south=xp.repeat(om_bc_south[:,None], repeats=nk, axis=1)
247
248 xwall_n=0.5*(x2d[0:-1,-1]+x2d[1:,-1])
249 ywall_n=0.5*(y2d[0:-1,-1]+y2d[1:,-1])
250 dist2_n=(yp2d[:,-1]-ywall_n)**2+(xp2d[:,-1]-xwall_n)**2
251 om_bc_north=6*viscos/c_omega_2/dist2_n
252
253# make it 2D
254 om_bc_north=xp.repeat(om_bc_north[:,None], repeats=nk, axis=1)
255
256 om_bc_west_type='d'
257 om_bc_east_type='n'
258 om_bc_south_type='d'
259 om_bc_north_type='d'
260 om_bc_low_type='d'
261 om_bc_high_type='d'
262
263 return
264
265
266
267def modify_init(u3d,v3d,w3d,k3d,om3d,eps3d,vis3d):
268
269 return u3d,v3d,w3d,k3d,om3d,eps3d,vis3d,dist3d
270
272 return u_bc_west,v_bc_west,w_bc_west,k_bc_west,eps_bc_west,om_bc_west,u3d_face_w,convw
273
274def modify_conv(convw,convs,convl):
275
276 return convw,convs,convl
277
278def modify_u(su3d,sp3d):
279
280 su3d= su3d+vol
281
282 if iter == 0:
283 if itstep == 0:
284 print('file opened')
285 xp.savetxt('u-time-history.dat',\
286 xp.c_[itstep,k3d[1,5,0],k3d[1,10,0],k3d[1,20,0],k3d[1,30,0],k3d[1,40,0],\
287 k3d[1,50,0],k3d[1,60,0]])
288 else:
289 print('file printed')
290 with open('u-time-history.dat','ab') as f:
291 xp.savetxt(f,\
292 xp.c_[itstep,k3d[1,5,0],k3d[1,10,0],k3d[1,20,0],k3d[1,30,0],k3d[1,40,0],\
293 k3d[1,50,0],k3d[1,60,0]])
294
295 return su3d,sp3d
296
297
298def modify_v(su3d,sp3d):
299
300 return su3d,sp3d
301
302
303def modify_w(su3d,sp3d):
304
305 return su3d,sp3d
306
307
308def modify_k(su3d,sp3d,gen):
309
310 comm_term=xp.zeros((ni,nj,nk))
311
312 return su3d,sp3d,comm_term
313
314
315def modify_eps(su3d,sp3d):
316
317 return su3d,sp3d
318
319def modify_om(su3d,sp3d,comm_term):
320
321 return su3d,sp3d
322
323def modify_outlet(convw):
324
325 return convw,u_bc_east
326
327def fix_k():
328
329 return aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d
330
332
333 return aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d
334
336
337 return aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d
338
339
340def modify_vis(vis3d):
341
342 return vis3d
343
344def STG(U0,viscos,itstep,deltat,xp2d_synt,yp2d_synt,zp2d_synt,dw2d,dx2d,dy2d,dz2d,r11,r12,r13,r22,r23,r33,k_rans,om_rans):
345#!!! number of modes = nmodes
346#!!! kinetic viscosity = visc
347
348 global N,kn_wave,le_max,kx,ky,kz,psi,qnNorm,sxio,syio,szio,a11,a21,a22,a31,a32,a33
349
350 nj = int(xp.size(yp2d_synt,0))
351 nk = int(xp.size(yp2d_synt,1))
352
353 if itstep == 0:
354
355 uvmean_check=0
356 uv_synt_mean=0
357# compute length scales le, lcut, lnu, lt
358 eps = 0.09*om_rans*k_rans
359
360 lt = k_rans**0.5/(0.09*om_rans)
361 le = xp.minimum(2.0*dw2d,3.0*lt)
362 lnu = viscos**0.75/eps**0.25
363 le_max = xp.amax(le)
364
365 tmp = xp.maximum(dy2d,dz2d)
366 hmax = xp.maximum(tmp,dx2d)
367 tmpNew = xp.maximum(tmp,0.3*hmax) + 0.1*dw2d
368 lcut = 2.0*xp.minimum(tmpNew,hmax)
369# compute ke, kcut, knu
370 ke = 2.0*xp.pi/le
371 kcut = 2.0*xp.pi/lcut
372 knu = 2*xp.pi/lnu
373
374 ke_min = 2*xp.pi/le_max
375 k_min_STG = 0.5*ke_min
376 k_max = 1.5*xp.amax(kcut)
377 alpha = 0.01
378
379 N = int(xp.ceil(xp.log(k_max/k_min_STG)/xp.log(1+alpha) + 1))
380 n = xp.linspace(1,N,N)
381 kn = k_min_STG*(1+alpha)**(n-1)
382 dkn = xp.zeros(N)
383 dkn[0] = 0.5*(kn[1]-kn[0])
384 dkn[1:-1] = 0.5*(kn[2:]-kn[0:-2])
385 dkn[-1] = 0.5*(kn[-1]-kn[-2])
386
387# fcut =
388
389# create a seed from time
390 xp.random.seed()
391 xp.random.seed(2)
392
393# zero all arrays to zero
394# wnr=xp.zeros(nmodes+2)
395# fi=xp.zeros(nmodes+2)
396 teta=xp.zeros(N)
397 psi=xp.zeros(N)
398# wnr=xp.zeros(nmodes+2)
399 kxio=xp.zeros((nj,nk,N))
400 kyio=xp.zeros((nj,nk,N))
401 kzio=xp.zeros((nj,nk,N))
402 sxio=xp.zeros((nj,nk,N))
403 syio=xp.zeros((nj,nk,N))
404 szio=xp.zeros((nj,nk,N))
405# yp2d_wave=xp.zeros((nj,nk,nmodes+2))
406 zp2d_wave=xp.zeros((nj,nk,N))
407 u=xp.zeros((nj,nk))
408 v=xp.zeros((nj,nk))
409 w=xp.zeros((nj,nk))
410 # # compute random angles
411
412 fi = xp.random.uniform(0.,2.*math.pi,N)
413 psi = xp.random.uniform(0.,2.*math.pi,N)
414 alfa = xp.random.uniform(0.,2.*math.pi,N)
415 ang = xp.random.uniform(0.,1,N)
416 teta=xp.arccos(1.-ang/0.5)
417
418# wavenumber vector from random angles
419 kxio=xp.sin(teta)*xp.cos(fi)
420 kyio=xp.sin(teta)*xp.sin(fi)
421 kzio=xp.cos(teta)
422#
423# sigma (s=sigma) from random angles. sigma is the unit direction which gives the direction
424# of the synthetic velocity vector (u, v, w)
425 sxio=xp.cos(fi)*xp.cos(teta)*xp.cos(alfa)-xp.sin(fi)*xp.sin(alfa)
426 syio=xp.sin(fi)*xp.cos(teta)*xp.cos(alfa)+xp.cos(fi)*xp.sin(alfa)
427 szio=-xp.sin(teta)*xp.cos(alfa)
428
429 kcut_wave = xp.repeat(kcut[:,:,None], repeats=N, axis=2)
430 knu_wave = xp.repeat(knu[:,:,None], repeats=N, axis=2)
431 ke_wave = xp.repeat(ke[:,:,None], repeats=N, axis=2)
432 kn_wave = xp.repeat(kn[:,None],repeats=nj,axis=1)
433 kn_wave = xp.repeat(kn_wave[:,:,None],repeats=nk,axis=2)
434 kn_wave = xp.transpose(kn_wave,(1,2,0))
435
436 fnu = xp.exp(-(12.0*kn_wave/knu_wave)**2)
437 fcut = xp.exp(-(4.0*xp.maximum(kn_wave-0.9*kcut_wave,xp.zeros((nj,nk,N)))/kcut_wave)**3)
438
439 E = (kn_wave/ke_wave)**4/(1.0+2.4*(kn_wave/ke_wave)**2)**(17/6)*fnu*fcut
440
441 # print('ASD')
442 # jidx = 35
443 # kidx = 15
444
445 # Eplot = E[jidx,kidx,:]
446
447 # fig1,ax1 = plt.subplots()
448 # plt.subplots_adjust(left=0.20,bottom=0.20)
449 # plt.loglog(kn,Eplot,'b--',label="$x=0$")
450 # plt.loglog(kn,kn**(-5/3),'r-',label="$x=0$")
451 # plt.axis([1, xp.amax(kn), 0.000001, 1])
452
453 qn = E*dkn
454 qnSum = xp.sum(qn,axis=2)
455
456 qnNorm = qn/xp.repeat(qnSum[:,:,None], repeats=N, axis=2)
457
458 kx=kxio*kn
459 ky=kyio*kn
460 kz=kzio*kn
461
462 # Cholesky Decomposition
463 a11 = xp.sqrt(r11)
464 a21 = r12/(a11)
465 a22 = xp.sqrt(r22 - a21*a21)
466 a31 = r13/(a11)
467 a32 = (r23 - a21*a31)/(a22)
468 a33 = xp.sqrt(r33 - a31*a31 - a32*a32)
469
470 print('le_min : ' + str(xp.amin(le)))
471 print('lcut_min : ' + str(xp.amin(lcut)))
472 print('lnu_min : ' + str(xp.amin(lnu)))
473 print('k_min_STG : ' + str(k_min_STG))
474 print('k_max : ' + str(k_max))
475 print('NModes STG: ' + str(N))
476
477# #
478# #=========================================================================
479# #
480
481 xp2d_wave=xp.repeat(xp2d_synt[:,:,None], repeats=N, axis=2)
482 arg1=(2*xp.pi/(kn_wave*le_max))*(xp2d_wave-U0*deltat*itstep)*kx
483
484 yp2d_wave=xp.repeat(yp2d_synt[:,:,None], repeats=N, axis=2)
485 arg2=yp2d_wave*ky
486
487 zp2d_wave=xp.repeat(zp2d_synt[:,:,None], repeats=N, axis=2)
488 arg3=zp2d_wave*kz
489
490 arg=arg1+arg2+arg3+psi
491
492 tfunk=xp.cos(arg)
493
494# sum over all wavenumbers => synthetic velocity field
495 usynt= xp.sqrt(6.0)*xp.sum(xp.sqrt(qnNorm)*tfunk*sxio,axis=2)
496 vsynt= xp.sqrt(6.0)*xp.sum(xp.sqrt(qnNorm)*tfunk*syio,axis=2)
497 wsynt= xp.sqrt(6.0)*xp.sum(xp.sqrt(qnNorm)*tfunk*szio,axis=2)
498
499#
500 usynt_aniso = a11*usynt
501 vsynt_aniso = a21*usynt+a22*vsynt
502 wsynt_aniso = a31*usynt+a32*vsynt+a33*wsynt
503
504 return usynt_aniso,vsynt_aniso,wsynt_aniso
505
506import cupy as xp
507import time,random,sys
508#from scipy.signal import welch, hann
509import math
510
511
512def synt_fluct(nmodes,it,sli,yp,zp,uv_rans,visc,jmirror,dmin_synt):
513#=========================== chapter 1 ============================================
514
515#!!! number of modes = nmodes
516#!!! smallest wavenumber = dxmin
517#!!! ratio of ke and kmin (in wavenumber) = wew1fct
518#!!! turb. velocity scale = up
519#!!! diss. rate. = epsm
520#!!! kinetic viscosity = visc
521#!!! length scale = sli
522#!!! mirror vfluct at j > jmirror
523
524 global dxmin,amp,epsm,epsm,wnr1,wew1fct,xp,yp2d_synt,zp2d_synt,xp2d_synt,nj,up
525 global e,kxio,kyio,kyio,sxio,syio,syio,utn,tfunk,wnre,dkn,arg1,arg2,arg3,arg,fi,psi,teta,alfa,wnr,kx,ky,kz,wnrn,\
526 r11,r12,r13,r21,r22,r23,r31,r32,r33,a11,a22,a33,wnreta,uv_synt_mean,uvmean_check,a11i,a22i,a33i,\
527 xp2d_wave,yp2d_wave,zp2d_wave,uv_rans_non,uv_rans_max
528 global utn,tfunk,sx,e,rk,kxi,usynt_wave,usynt1,usynt,sy,vsynt,yp2d_org,usynt_aniso,scale_2
529
530 uv_rans=xp.abs(uv_rans)
531 if it == 0:
532
533 uvmean_check=0
534 uv_synt_mean=0
535# anisotropix fluctuations
536# R=xp.loadtxt('R.dat')
537 R=xp.genfromtxt("R.dat", dtype=None,comments="%")
538 r11=R[0,0]
539 r12=R[0,1]
540 r13=R[0,2]
541 r21=R[1,0]
542 r22=R[1,1]
543 r23=R[1,2]
544 r31=R[2,0]
545 r32=R[2,1]
546 r33=R[2,2]
547
548# A=xp.loadtxt('a.dat')
549 A=xp.genfromtxt("a.dat", dtype=None,comments="%")
550 a11=A[0]
551 a22=A[1]
552 a33=A[2]
553
554 amp=1.452762113
555 wew1fct=2
556
557# in log region: k/uvmax=3.3 => up=(3.3*uvmax)**0.5
558 if xp.all(uv_rans==1):
559 up=1
560 else:
561 up=(3.3*xp.max(uv_rans))**0.5
562
563 epsm=up**3/sli
564#
565#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
566# number of grid points in y, z
567 nj=len(yp)
568 nk=len(zp)
569
570# make it 2D
571 zp2d_synt=xp.repeat(zp[None,:], repeats=nj, axis=0)
572
573# make it 2D
574 yp2d_synt=xp.repeat(yp[:,None], repeats=nk, axis=1)
575
576
577 xpp=0.5
578# make it 2D
579 xp2d_synt=xp.ones((nj,nk))*xpp
580
581 yp2d_org=yp2d_synt
582# transform to principal coord. directions
583 xp2d_synt=r11*xpp+r21*yp2d_org+r31*zp2d_synt
584 yp2d_synt=r12*xpp+r22*yp2d_synt+r32*zp2d_synt
585 zp2d_synt=r13*xpp+r23*yp2d_synt+r33*zp2d_synt
586
587# search min grid step
588 dminy=xp.min(xp.diff(yp))
589 dminz=xp.min(xp.diff(zp))
590 dxmin=min(dminy,dminz)
591
592# don't let is be smaller than dmin_synt
593 dxmin=max(dxmin,dmin_synt)
594
595
596
597# create a seed from time
598 xp.random.seed()
599 xp.random.seed(2)
600
601# zero all arrays to zero
602 wnr=xp.zeros(nmodes+2)
603 fi=xp.zeros(nmodes+2)
604 teta=xp.zeros(nmodes+2)
605 psi=xp.zeros(nmodes+2)
606 wnr=xp.zeros(nmodes+2)
607 kxio=xp.zeros((nj,nk,nmodes+2))
608 kyio=xp.zeros((nj,nk,nmodes+2))
609 kzio=xp.zeros((nj,nk,nmodes+2))
610 sxio=xp.zeros((nj,nk,nmodes+2))
611 syio=xp.zeros((nj,nk,nmodes+2))
612 szio=xp.zeros((nj,nk,nmodes+2))
613# yp2d_wave=xp.zeros((nj,nk,nmodes+2))
614 zp2d_wave=xp.zeros((nj,nk,nmodes+2))
615 u=xp.zeros((nj,nk))
616 v=xp.zeros((nj,nk))
617 w=xp.zeros((nj,nk))
618#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
619
620# highest wave number
621 wnrn=2.*math.pi/dxmin
622#
623# k_e (related to peak energy wave number)
624 wnre=9.*math.pi*amp/(55.*sli)
625#
626# wavenumber used in the viscous expression (high wavenumbers) in the von Karman spectrum
627 wnreta=(epsm/visc**3)**0.25
628
629# smallest wavenumber
630 wnr1=wnre/wew1fct
631
632# wavenumber step
633 dkn=(wnrn-wnr1)/nmodes
634
635# wavenumers
636 wnr=xp.linspace(wnr1,wnrn,nmodes)
637
638# invert the eigenvalue matrix (anisotropic)
639 a11i=1/a11
640 a22i=1/a22
641 a33i=1/a33
642
643# make a non-dimensional uv_rans profile
644 uv_rans_max=xp.max(uv_rans)
645 uv_rans_non=(uv_rans/(uv_rans_max+1e-10))**0.5
646# make it 2D
647 uv_rans_non=xp.repeat(uv_rans_non[:,None], repeats=nk, axis=1)
648
649
650 print(f"\n{'sli: '} {sli}, {'visc: '}{visc:.2e}, {'nmodes: '}{nmodes}, {'dxmin: '}{dxmin:.3e}, {'dkn: '}{dkn:.3e}, {'dmin_synt: '}{dmin_synt:.3e}")
651
652 print(f"\n{'wnre: '} {wnre:.2e}, {'wnr1: '}{wnr1:.2e}, {'epsm: '}{epsm:.3e}, {'wnrn: '}{wnrn:.3e}")
653
654 print(f"\n{'eigenvalue 1, 2 and 3: '}{a11:.3e}, {a22:.3e}, {a33:.3e}")
655 print(f"\n{'eigenvector R11, R12 and R13: '}{r11:.3e}, {r12:.3e}, {r13:.3e}")
656 print(f"\n{'eigenvector R21, R22 and R23: '}{r21:.3e}, {r22:.3e}, {r23:.3e}")
657 print(f"\n{'eigenvector R31, R32 and R33: '}{r31:.3e}, {r32:.3e}, {r33:.3e}\n")
658
659
660
661#
662#=========================== chapter 2 ============================================
663#
664
665# compute random angles
666 fi = xp.random.uniform(0.,2.*math.pi,nmodes)
667 psi = xp.random.uniform(0.,2.*math.pi,nmodes)
668 alfa = xp.random.uniform(0.,2.*math.pi,nmodes)
669 ang = xp.random.uniform(0.,1,nmodes)
670 teta=xp.arccos(1.-ang/0.5)
671
672 print('time step no,',it)
673
674
675# wavenumber vector from random angles
676 kxio=xp.sin(teta)*xp.cos(fi)
677 kyio=xp.sin(teta)*xp.sin(fi)
678 kzio=xp.cos(teta)
679#
680# sigma (s=sigma) from random angles. sigma is the unit direction which gives the direction
681# of the synthetic velocity vector (u, v, w)
682 sxio=xp.cos(fi)*xp.cos(teta)*xp.cos(alfa)-xp.sin(fi)*xp.sin(alfa)
683 syio=xp.sin(fi)*xp.cos(teta)*xp.cos(alfa)+xp.cos(fi)*xp.sin(alfa)
684 szio=-xp.sin(teta)*xp.cos(alfa)
685
686#
687#=========================== chapter 3 ============================================
688#
689# loop over all wavenumbers
690 kxi=r11*kxio+r21*kyio+r31*kzio
691 kyi=r12*kxio+r22*kyio+r32*kzio
692 kzi=r13*kxio+r23*kyio+r33*kzio
693
694 sxi=r11*sxio+r21*syio+r31*szio
695 syi=r12*sxio+r22*syio+r32*szio
696 szi=r13*sxio+r23*syio+r33*szio
697
698 sx=a11**0.5*sxi
699 sy=a22**0.5*syi
700 sz=a33**0.5*szi
701
702 kx=kxi*wnr*a11i**0.5
703 ky=kyi*wnr*a22i**0.5
704 kz=kzi*wnr*a33i**0.5
705 rk=xp.sqrt(kx**2+ky**2+kz**2)
706
707
708
709 xp2d_wave=xp.repeat(xp2d_synt[:,:,None], repeats=nmodes, axis=2)
710 arg1=xp2d_wave*kx
711
712
713 yp2d_wave=xp.repeat(yp2d_synt[:,:,None], repeats=nmodes, axis=2)
714 arg2=yp2d_wave*ky
715
716 zp2d_wave=xp.repeat(zp2d_synt[:,:,None], repeats=nmodes, axis=2)
717 arg3=zp2d_wave*kz
718
719 arg=arg1+arg2+arg3+psi
720
721 tfunk=xp.cos(arg)
722
723# von Karman spectrum
724 e=amp/wnre*(wnr/wnre)**4/((1.+(wnr/wnre)**2)**(17./6.))*xp.exp(-2*(wnr/wnreta)**2)
725
726# include only wavenumber for which rk < wnrn
727 e=xp.where(rk < wnrn,e,0)
728
729 utn=xp.sqrt(e*up**2*dkn)
730
731# sum over all wavenumbers => synthetic velocity field
732 usynt=xp.sum(2.*utn*tfunk*sx,axis=2)
733 vsynt=xp.sum(2.*utn*tfunk*sy,axis=2)
734 wsynt=xp.sum(2.*utn*tfunk*sz,axis=2)
735
736# transform back to x-y-z => anjsotropic fluct
737 usynt_aniso=r11*usynt+r12*vsynt+r13*wsynt
738 vsynt_aniso=r21*usynt+r22*vsynt+r23*wsynt
739 wsynt_aniso=r31*usynt+r32*vsynt+r33*wsynt
740
741# mean shear stress (must be computed before mirroring)
742 uv=xp.mean(usynt_aniso*vsynt_aniso)
743
744# mirror vfluct
745 if jmirror > 0:
746 vsynt_aniso[jmirror:,:]=-vsynt_aniso[jmirror:,:]
747
748# sum over timesteps
749 uv_synt_mean=uv_synt_mean+uv
750
751# compute average
752 uvmean_time=xp.abs(uv_synt_mean)/(it+1)
753 print('uvmean_time',uvmean_time)
754
755 scale_2=(uv_rans_max/uvmean_time)**0.5*uv_rans_non
756
757# if uv_rans=1: don't scale
758# scale_2=xp.where(uv_rans==1,1,(uv_rans_max/uvmean_time)**0.5*uv_rans_non)
759
760# scale all fluctuations with uv_rans
761 usynt_aniso=usynt_aniso*scale_2
762 vsynt_aniso=vsynt_aniso*scale_2
763 wsynt_aniso=wsynt_aniso*scale_2
764
765
766# compute mean of synt fluct
767 uvmean_check=uvmean_check+xp.mean(usynt_aniso*vsynt_aniso,axis=1)
768
769# peak of uv_rans
770 j=xp.where(uv_rans == xp.amax(uv_rans))
771
772# check peak
773 print('synt: uvmean',xp.abs(uvmean_check[j])/(it+1),'uv_rans_max=',xp.max(xp.abs(uv_rans)),'at j=',j)
774
775 return usynt_aniso,vsynt_aniso,wsynt_aniso
776if gpu:
777 from cupyx.scipy import sparse
778 import cupy as xp
779 from cupyx.scipy.sparse import spdiags,linalg,eye
780else:
781 from scipy import sparse
782 import numpy as xp
783 from scipy.sparse import spdiags,linalg,eye
784import sys
785import time
786import pyamg
787import pyamgx
788
789import socket
790
791def init():
792 print('hostname: ',socket.gethostname())
793
794# distance to nearest wall
795 ywall_s=0.5*(y2d[0:-1,0]+y2d[1:,0])
796 dist_s=yp2d-ywall_s[:,None]
797 ywall_n=0.5*(y2d[0:-1,-1]+y2d[1:,-1])
798 dist_n=ywall_n[:,None] -yp2d
799 dist=xp.minimum(dist_s,dist_n)
800 dist3d=xp.repeat(dist[:,:,None], repeats=nk, axis=2)
801
802# west face coordinate
803 xw=0.5*(x2d[0:-1,0:-1]+x2d[0:-1,1:])
804 yw=0.5*(y2d[0:-1,0:-1]+y2d[0:-1,1:])
805
806 del1x=((xw-xp2d)**2+(yw-yp2d)**2)**0.5
807 del2x=((xw-xp.roll(xp2d,1,axis=0))**2+(yw-xp.roll(yp2d,1,axis=0))**2)**0.5
808 fx=del2x/(del1x+del2x)
809 fx = xp.dstack([fx]*nk)
810
811 if cyclic_x:
812 fx[0,:,:]=0.5
813
814# south face coordinate
815 xs=0.5*(x2d[0:-1,0:-1]+x2d[1:,0:-1])
816 ys=0.5*(y2d[0:-1,0:-1]+y2d[1:,0:-1])
817
818 del1y=((xs-xp2d)**2+(ys-yp2d)**2)**0.5
819 del2y=((xs-xp.roll(xp2d,1,axis=1))**2+(ys-xp.roll(yp2d,1,axis=1))**2)**0.5
820 fy=del2y/(del1y+del2y)
821 fy = xp.dstack([fy]*nk)
822
823# low face coordinate
824 zl=z[0:-1]
825
826 del1z=zp-zl
827 del2z=xp.abs(zl-xp.roll(zp,1))
828 fz=del2z/(del1z+del2z)
829 fz=xp.repeat(fz[None,:], repeats=nj, axis=0)
830 fz=xp.repeat(fz[None,:,:], repeats=ni, axis=0)
831
832 if cyclic_z:
833 fz[:,:,0]=0.5
834
835 areawy=xp.diff(x2d,axis=1)
836 areawx=-xp.diff(y2d,axis=1)
837
838# make them 3d
839 areawx= xp.dstack([areawx]*nk)
840 areawy= xp.dstack([areawy]*nk)
841 areawy=areawy*dz3d[:,1:,:]
842 areawx=areawx*dz3d[:,1:,:]
843
844 areasy=-xp.diff(x2d,axis=0)
845 areasx=xp.diff(y2d,axis=0)
846# make them 3d
847 areasx= xp.dstack([areasx]*nk)
848 areasy= xp.dstack([areasy]*nk)
849 areasy=areasy*dz3d[1:,:,:]
850 areasx=areasx*dz3d[1:,:,:]
851
852 areaw=(areawx**2+areawy**2)**0.5
853 areas=(areasx**2+areasy**2)**0.5
854
855# volume approaximated as the vector product of two triangles for cells
856 ax=xp.diff(x2d,axis=1)
857 ay=xp.diff(y2d,axis=1)
858 bx=xp.diff(x2d,axis=0)
859 by=xp.diff(y2d,axis=0)
860
861 areaz_1=0.5*xp.absolute(ax[0:-1,:]*by[:,0:-1]-ay[0:-1,:]*bx[:,0:-1])
862
863 ax=xp.diff(x2d,axis=1)
864 ay=xp.diff(y2d,axis=1)
865 areaz_2=0.5*xp.absolute(ax[1:,:]*by[:,0:-1]-ay[1:,:]*bx[:,0:-1])
866
867 areaz=areaz_1+areaz_2
868# make it 3d
869 areal= xp.dstack([areaz]*(nk+1))
870 vol=areal[:,:,1:]*dz3d[0:-1,0:-1,:]
871
872# coeff at south wall (without viscosity)
873 as_bound=areas[:,0,:]**2/(0.5*vol[:,0,:])
874
875# coeff at north wall (without viscosity)
876 an_bound=areas[:,-1,:]**2/(0.5*vol[:,-1,:])
877
878# coeff at west wall (without viscosity)
879 aw_bound=areaw[0,:,:]**2/(0.5*vol[0,:,:])
880 if cyclic_x:
881 aw_bound=areaw[0,:,:]**2/(0.5*(vol[0,:,:]+vol[-1,:,:]))
882
883# coeff at east wall (without viscosity) N.B: which cyclic_x
884# this is never used
885 ae_bound=areaw[-1,:,:]**2/(0.5*vol[-1,:,:])
886
887# make it 2d
888 al_bound=areal[:,:,0]**2/(0.5*vol[:,:,0])
889 if cyclic_z:
890 al_bound=areal[:,:,0]**2/(0.5*(vol[:,:,0]+vol[:,:,-1]))
891
892 ah_bound=areal[:,:,-1]**2/(0.5*vol[:,:,-1])
893
894 return areaw,areawx,areawy,areas,areasx,areasy,areal,vol,fx,fy,fz,aw_bound,ae_bound,as_bound,an_bound,al_bound,ah_bound,dist3d
895
897
898 print('////////////////// Start of input data ////////////////// \n\n\n')
899
900 print('\n\n########### section 0 section 0 choice of CPU or GPU ##########')
901 print(f"{'GPU: ':<29} {gpu}")
902
903 print('\n\n########### section 1 choice of differencing scheme ###########')
904 print(f"{'scheme: ':<29} {scheme}")
905 print(f"{'scheme_turb: ':<29} {scheme_turb}")
906 print(f"{'acrank: ':<29} {acrank}")
907 print(f"{'acrank_conv: ':<29} {acrank_conv}")
908 print(f"{'acrank_conv_keps: ':<29} {acrank_conv_keps}")
909 print(f"{'acrank_conv_kom: ':<29} {acrank_conv_kom}")
910 if scheme == 'm':
911 print(f"{'blend: ':<29} {blend}")
912
913
914 print('\n\n########### section 2 turbulence models ###########')
915
916 print(f"{'cmu: ':<29} {cmu}")
917 print(f"{'pans: ':<29} {pans}")
918 print(f"{'keps: ':<29} {keps}")
919 print(f"{'kom_des: ':<29} {kom_des}")
920 print(f"{'keps_des: ':<29} {keps_des}")
921 print(f"{'k_eq_les: ':<29} {k_eq_les}")
922 if jl0 < 0:
923 print(f"{'jl0: ':<29} {jl0}")
924 print(f"{'kom: ':<29} {kom}")
925 print(f"{'sst: ':<29} {sst}")
926 print(f"{'smag: ':<29} {smag}")
927 print(f"{'wale: ':<29} {wale}")
928 if sst:
929 print(f"{'prand_k_sst_1: ':<29} {prand_k_sst_1}")
930 print(f"{'prand_k_sst_2: ':<29} {prand_k_sst_2}")
931 print(f"{'prand_omega_sst_1: ':<29} {prand_omega_sst_1}")
932 print(f"{'prand_omega_sst_2: ':<29} {prand_omega_sst_2}")
933 print(f"{'cdes: ':<29} {cdes}")
934 if k_eq_les:
935 print(f"{'c_eps: ':<29} {c_eps}")
936 if pans:
937 print(f"{'fkmin_limit: ':<29} {fkmin_limit}")
938 if sst :
939 print(f"{'c_omega_1_sst_1: ':<29} {c_omega_1_sst_1:.3f}")
940 print(f"{'c_omega_1_sst_2: ':<29} {c_omega_1_sst_2:.3f}")
941 print(f"{'c_omega_2_sst_1: ':<29} {c_omega_2_sst_1}")
942 print(f"{'c_omega_2_sst_2: ':<29} {c_omega_2_sst_2}")
943 if keps or pans or keps_des:
944 print(f"{'c_eps_1: ':<29} {c_eps_1}")
945 print(f"{'c_eps_2: ':<29} {c_eps_2}")
946 print(f"{'prand_k: ':<29} {prand_k}")
947 print(f"{'prand_eps: ':<29} {prand_eps}")
948 if kom or kom_des:
949 print(f"{'c_omega_1: ':<29} {c_omega_1:.3f}")
950 print(f"{'c_omega_2: ':<29} {c_omega_2}")
951 print(f"{'prand_k: ':<29} {prand_k}")
952 print(f"{'prand_omega: ':<29} {prand_omega}")
953
954 if keps_des or kom_des:
955 print(f"{'cdes: ':<29} {cdes}")
956
957
958 print('\n\n########### section 3 restart/save ###########')
959 print(f"{'restart: ':<29} {restart}")
960 print(f"{'save: ':<29} {save}")
961
962
963 print('\n\n########### section 4 fluid properties ###########')
964 print(f"{'viscos: ':<29} {viscos:.2e}")
965
966 print('\n\n########### section 5 relaxation factors ###########')
967 print(f"{'urfvis: ':<29} {urfvis}")
968
969
970 print('\n\n########### section 6 number of iteration and convergence criterira ###########')
971 print(f"{'sormax: ':<29} {sormax}")
972 print(f"{'norm_order: ':<29} {norm_order}")
973 print(f"{'maxit: ':<29} {maxit}")
974 print(f"{'min_iter: ':<29} {min_iter}")
975 print(f"{'solver_vel: ':<29} {solver_vel}")
976 print(f"{'solver_p: ':<29} {solver_p}")
977 print(f"{'coeff_v: ':<29} {coeff_v}")
978 print(f"{'coeff_w: ':<29} {coeff_w}")
979 print(f"{'embedded: ':<29} {embedded}")
980 if embedded:
981 print(f"{'x_embed: ':<29} {x_embed}")
982 print(f"{'solver_turb: ':<29} {solver_turb}")
983 print(f"{'amg_relax: ':<29} {amg_relax}")
984 print(f"{'amg_cycle: ':<29} {amg_cycle}")
985 if solver_vel == 'pyamg' or solver_turb == 'pyamg':
986 if amg_relax_phi != 'default':
987 print(f"{'amg_relax_phi: ':<29} {amg_relax_phi}")
988 print(f"{'amg_cycle_phi: ':<29} {amg_cycle_phi}")
989 print(f"{'nsweep_vel: ':<29} {nsweep_vel}")
990 print(f"{'nsweep_keps: ':<29} {nsweep_keps}")
991 print(f"{'nsweep_kom: ':<29} {nsweep_kom}")
992 if gpu:
993 print(f"{'convergence_limit_gpu: ':<29} {convergence_limit_gpu}")
994 else:
995 print(f"{'convergence_limit_u: ':<29} {convergence_limit_u}")
996 print(f"{'convergence_limit_v: ':<29} {convergence_limit_v}")
997 print(f"{'convergence_limit_w: ':<29} {convergence_limit_w}")
998 print(f"{'convergence_limit_p: ':<29} {convergence_limit_p}")
999 if keps or pans or sst or kom or kom_des or keps_des:
1000 print(f"{'convergence_limit_k: ':<29} {convergence_limit_k}")
1001 if keps or pans or keps_des:
1002 print(f"{'convergence_limit_eps: ':<29} {convergence_limit_eps}")
1003 if sst or kom or kom_des:
1004 print(f"{'convergence_limit_om: ':<29} {convergence_limit_om}")
1005
1006 print('\n\n########### section 7 all variables are printed during the iteration at node ###########')
1007 print(f"{'imon: ':<29} {imon}")
1008 print(f"{'jmon: ':<29} {jmon}")
1009 print(f"{'kmon: ':<29} {kmon}")
1010
1011
1012 print('\n\n########### section 8 time-averaging ###########')
1013 print(f"{'ntstep: ':<29} {ntstep}")
1014 print(f"{'dt[0]: ':<29} {dt[0]:.2e}")
1015 print(f"{'itstep_start: ':<29} {itstep_start}")
1016 print(f"{'itstep_save: ':<29} {itstep_save}")
1017 print(f"{'itstep_stats: ':<29} {itstep_stats}")
1018
1019
1020 print('\n\n########### section 9 residual scaling parameters ###########')
1021 print(f"{'resnorm_p: ':<29} {resnorm_p:.1f}")
1022 print(f"{'resnorm_vel: ':<29} {resnorm_vel:.1f}")
1023
1024
1025 print('\n\n########### Section 10 grid and boundary conditions ###########')
1026 print(f"{'ni: ':<29} {ni}")
1027 print(f"{'nj: ':<29} {nj}")
1028 print(f"{'nk: ':<29} {nk}")
1029 if xp.isscalar(dz):
1030 print(f"{'dz: ':<29} {dz:.2e}")
1031 print('\n')
1032 print(f"{'cyclic_x: ':<29} {cyclic_x}")
1033 print(f"{'cyclic_z: ':<29} {cyclic_z}")
1034 print('\n')
1035 if not cyclic_x:
1036 print(f"{'L_t_synt: ':<29} {L_t_synt}")
1037 print(f"{'nmodes_synt: ':<29} {nmodes_synt}")
1038 print(f"{'dmin_synt: ':<29} {dmin_synt}")
1039 print(f"{'jmirror_synt: ':<29} {jmirror_synt}")
1040 print('\n')
1041
1042 if i_block_start !=0 or i_block_end !=0 or j_block_start !=0 or i_block_end !=0:
1043 print('------blockage')
1044 print(f"{' ':<5}{'i_block_start: ':<29} {i_block_start}")
1045 print(f"{' ':<5}{'i_block_end: ':<29} {i_block_end}")
1046 print(f"{' ':<5}{'j_block_start: ':<29} {j_block_start}")
1047 print(f"{' ':<5}{'j_block_end: ':<29} {j_block_end}")
1048 print('\n')
1049
1050 print('------boundary conditions for u')
1051 if not cyclic_x:
1052 print(f"{' ':<5}{'u_bc_west_type: ':<29} {u_bc_west_type}")
1053 print(f"{' ':<5}{'u_bc_east_type: ':<29} {u_bc_east_type}")
1054 if u_bc_west_type == 'd':
1055 print(f"{' ':<5}{'u_bc_west[0,0]: ':<29} {u_bc_west[0,0]}")
1056 if u_bc_east_type == 'd':
1057 print(f"{' ':<5}{'u_bc_east[0,0]: ':<29} {u_bc_east[0,0]}")
1058
1059
1060 print(f"{' ':<5}{'u_bc_south_type: ':<29} {u_bc_south_type}")
1061 print(f"{' ':<5}{'u_bc_north_type: ':<29} {u_bc_north_type}")
1062
1063 if u_bc_south_type == 'd':
1064 print(f"{' ':<5}{'u_bc_south[0,0]: ':<29} {u_bc_south[0,0]}")
1065 if u_bc_north_type == 'd':
1066 print(f"{' ':<5}{'u_bc_north[0,0]: ':<29} {u_bc_north[0,0]}")
1067
1068 if not cyclic_z:
1069 print(f"{' ':<5}{'u_bc_low_type: ':<29} {u_bc_low_type}")
1070 print(f"{' ':<5}{'u_bc_high_type: ':<29} {u_bc_high_type}")
1071 if u_bc_low_type == 'd':
1072 print(f"{' ':<5}{'u_bc_lo[0,0]: ':<29} {u_bc_low[0,0]}")
1073 if u_bc_high_type == 'd':
1074 print(f"{' ':<5}{'u_bc_hig[0,0]: ':<29} {u_bc_high[0,0]}")
1075
1076 print('------boundary conditions for v')
1077 if not cyclic_x:
1078 print(f"{' ':<5}{'v_bc_west_type: ':<29} {v_bc_west_type}")
1079 print(f"{' ':<5}{'v_bc_east_type: ':<29} {v_bc_east_type}")
1080 if v_bc_west_type == 'd':
1081 print(f"{' ':<5}{'v_bc_west[0,0]: ':<29} {v_bc_west[0,0]}")
1082 if v_bc_east_type == 'd':
1083 print(f"{' ':<5}{'v_bc_east[0,0]: ':<29} {v_bc_east[0,0]}")
1084
1085
1086 print(f"{' ':<5}{'v_bc_south_type: ':<29} {v_bc_south_type}")
1087 print(f"{' ':<5}{'v_bc_north_type: ':<29} {v_bc_north_type}")
1088
1089 if v_bc_south_type == 'd':
1090 print(f"{' ':<5}{'v_bc_south[0,0]: ':<29} {v_bc_south[0,0]}")
1091 if v_bc_north_type == 'd':
1092 print(f"{' ':<5}{'v_bc_north[0,0]: ':<29} {v_bc_north[0,0]}")
1093
1094 if not cyclic_z:
1095 print(f"{' ':<5}{'v_bc_low_type: ':<29} {v_bc_low_type}")
1096 print(f"{' ':<5}{'v_bc_high_type: ':<29} {v_bc_high_type}")
1097 if v_bc_low_type == 'd':
1098 print(f"{' ':<5}{'v_bc_lo[0,0]: ':<29} {v_bc_low[0,0]}")
1099 if v_bc_high_type == 'd':
1100 print(f"{' ':<5}{'v_bc_hig[0,0]: ':<29} {v_bc_high[0,0]}")
1101
1102 print('------boundary conditions for w')
1103 if not cyclic_x:
1104 print(f"{' ':<5}{'w_bc_west_type: ':<29} {w_bc_west_type}")
1105 print(f"{' ':<5}{'w_bc_east_type: ':<29} {w_bc_east_type}")
1106 if w_bc_west_type == 'd':
1107 print(f"{' ':<5}{'w_bc_west[0,0]: ':<29} {w_bc_west[0,0]}")
1108 if w_bc_east_type == 'd':
1109 print(f"{' ':<5}{'w_bc_east[0,0]: ':<29} {w_bc_east[0,0]}")
1110
1111
1112 print(f"{' ':<5}{'w_bc_south_type: ':<29} {w_bc_south_type}")
1113 print(f"{' ':<5}{'w_bc_north_type: ':<29} {w_bc_north_type}")
1114
1115 if w_bc_south_type == 'd':
1116 print(f"{' ':<5}{'w_bc_south[0,0]: ':<29} {w_bc_south[0,0]}")
1117 if w_bc_north_type == 'd':
1118 print(f"{' ':<5}{'w_bc_north[0,0]: ':<29} {w_bc_north[0,0]}")
1119
1120 if not cyclic_z:
1121 print(f"{' ':<5}{'w_bc_low_type: ':<29} {w_bc_low_type}")
1122 print(f"{' ':<5}{'w_bc_high_type: ':<29} {w_bc_high_type}")
1123 if w_bc_low_type == 'd':
1124 print(f"{' ':<5}{'w_bc_low[0,0]: ':<29} {w_bc_low[0,0]}")
1125 if w_bc_high_type == 'd':
1126 print(f"{' ':<5}{'w_bc_hig[0,0]: ':<29} {w_bc_high[0,0]}")
1127
1128 print('------boundary conditions for p')
1129 if not cyclic_x:
1130 print(f"{' ':<5}{'p_bc_west_type: ':<29} {p_bc_west_type}")
1131 print(f"{' ':<5}{'p_bc_east_type: ':<29} {p_bc_east_type}")
1132 if p_bc_west_type == 'd':
1133 print(f"{' ':<5}{'p_bc_west[0,0]: ':<29} {p_bc_west[0,0]}")
1134 if p_bc_east_type == 'd':
1135 print(f"{' ':<5}{'p_bc_east[0,0]: ':<29} {p_bc_east[0,0]}")
1136
1137
1138 print(f"{' ':<5}{'p_bc_south_type: ':<29} {p_bc_south_type}")
1139 print(f"{' ':<5}{'p_bc_north_type: ':<29} {p_bc_north_type}")
1140
1141 if p_bc_south_type == 'd':
1142 print(f"{' ':<5}{'p_bc_south[0,0]: ':<29} {p_bc_south[0,0]}")
1143 if p_bc_north_type == 'd':
1144 print(f"{' ':<5}{'p_bc_north[0,0]: ':<29} {p_bc_north[0,0]}")
1145
1146 if not cyclic_z:
1147 print(f"{' ':<5}{'p_bc_low_type: ':<29} {p_bc_low_type}")
1148 print(f"{' ':<5}{'p_bc_high_type: ':<29} {p_bc_high_type}")
1149 if p_bc_low_type == 'd':
1150 print(f"{' ':<5}{'p_bc_lo[0,0]: ':<29} {p_bc_low[0,0]}")
1151 if p_bc_high_type == 'd':
1152 print(f"{' ':<5}{'p_bc_hig[0,0]: ':<29} {p_bc_high[0,0]}")
1153
1154 if sst or kom or kom_des or keps or pans or keps_des or k_eq_les:
1155 print('------boundary conditions for k')
1156 if not cyclic_x:
1157 print(f"{' ':<5}{'k_bc_west_type: ':<29} {k_bc_west_type}")
1158 print(f"{' ':<5}{'k_bc_east_type: ':<29} {k_bc_east_type}")
1159 if k_bc_west_type == 'd':
1160 print(f"{' ':<5}{'k_bc_west[0,0]: ':<29} {k_bc_west[0,0]}")
1161 if k_bc_east_type == 'd':
1162 print(f"{' ':<5}{'k_bc_east[0,0]: ':<29} {k_bc_east[0,0]}")
1163
1164
1165 print(f"{' ':<5}{'k_bc_south_type: ':<29} {k_bc_south_type}")
1166 print(f"{' ':<5}{'k_bc_north_type: ':<29} {k_bc_north_type}")
1167
1168 if k_bc_south_type == 'd':
1169 print(f"{' ':<5}{'k_bc_south[0,0]: ':<29} {k_bc_south[0,0]}")
1170 if k_bc_north_type == 'd':
1171 print(f"{' ':<5}{'k_bc_north[0,0]: ':<29} {k_bc_north[0,0]}")
1172
1173 if not cyclic_z:
1174 print(f"{' ':<5}{'k_bc_low_type: ':<29} {k_bc_low_type}")
1175 print(f"{' ':<5}{'k_bc_high_type: ':<29} {k_bc_high_type}")
1176 if k_bc_low_type == 'd':
1177 print(f"{' ':<5}{'k_bc_lo[0,0]: ':<29} {k_bc_low[0,0]}")
1178 if k_bc_high_type == 'd':
1179 print(f"{' ':<5}{'k_bc_hig[0,0]: ':<29} {k_bc_high[0,0]}")
1180
1181
1182 if keps or pans or keps_des:
1183 print('------boundary conditions for eps')
1184 if not cyclic_x:
1185 print(f"{' ':<5}{'eps_bc_west_type: ':<29} {eps_bc_west_type}")
1186 print(f"{' ':<5}{'eps_bc_east_type: ':<29} {eps_bc_east_type}")
1187 if eps_bc_west_type == 'd':
1188 print(f"{' ':<5}{'eps_bc_west[0,0]: ':<29} {eps_bc_west[0,0]}")
1189 if eps_bc_east_type == 'd':
1190 print(f"{' ':<5}{'eps_bc_east[0,0]: ':<29} {eps_bc_east[0,0]}")
1191
1192
1193 print(f"{' ':<5}{'eps_bc_south_type: ':<29} {eps_bc_south_type}")
1194 print(f"{' ':<5}{'eps_bc_north_type: ':<29} {eps_bc_north_type}")
1195
1196 if eps_bc_south_type == 'd':
1197 print(f"{' ':<5}{'eps_bc_south[0,0]: ':<29} {eps_bc_south[0,0]}")
1198 if eps_bc_north_type == 'd':
1199 print(f"{' ':<5}{'eps_bc_north[0,0]: ':<29} {eps_bc_north[0,0]}")
1200
1201 if not cyclic_z:
1202 print(f"{' ':<5}{'eps_bc_low_type: ':<29} {eps_bc_low_type}")
1203 print(f"{' ':<5}{'eps_bc_high_type: ':<29} {eps_bc_high_type}")
1204 if eps_bc_low_type == 'd':
1205 print(f"{' ':<5}{'eps_bc_lo[0,0]: ':<29} {eps_bc_low[0,0]}")
1206 if eps_bc_high_type == 'd':
1207 print(f"{' ':<5}{'eps_bc_hig[0,0]: ':<29} {eps_bc_high[0,0]}")
1208
1209 if sst or kom or kom_des:
1210 print('------boundary conditions for omega')
1211 if not cyclic_x:
1212 print(f"{' ':<5}{'om_bc_west_type: ':<29} {om_bc_west_type}")
1213 print(f"{' ':<5}{'om_bc_east_type: ':<29} {om_bc_east_type}")
1214 if om_bc_west_type == 'd':
1215 print(f"{' ':<5}{'om_bc_west[0,0]: ':<29} {om_bc_west[0,0]:.1f}")
1216 if om_bc_east_type == 'd':
1217 print(f"{' ':<5}{'om_bc_east[0,0]: ':<29} {om_bc_east[0,0]:.1f}")
1218
1219
1220 print(f"{' ':<5}{'om_bc_south_type: ':<29} {om_bc_south_type}")
1221 print(f"{' ':<5}{'om_bc_north_type: ':<29} {om_bc_north_type}")
1222
1223 if om_bc_south_type == 'd':
1224 print(f"{' ':<5}{'om_bc_south[0,0]: ':<29} {om_bc_south[0,0]:.1f}")
1225 if om_bc_north_type == 'd':
1226 print(f"{' ':<5}{'om_bc_north[0,0]: ':<29} {om_bc_north[0,0]:.1f}")
1227
1228 if not cyclic_z:
1229 print(f"{' ':<5}{'om_bc_low_type: ':<29} {om_bc_low_type}")
1230 print(f"{' ':<5}{'om_bc_high_type: ':<29} {om_bc_high_type}")
1231 if om_bc_low_type == 'd':
1232 print(f"{' ':<5}{'om_bc_low[0,0]: ':<29} {om_bc_low[0,0]}")
1233 if om_bc_high_type == 'd':
1234 print(f"{' ':<5}{'om_bc_high[0,0]: ':<29} {om_bc_high[0,0]}")
1235
1236
1237 print('\n\n\n ////////////////// End of input data //////////////////\n\n\n')
1238
1239
1240
1241 return
1242
1243def compute_face_phi(phi3d,phi_bc_west,phi_bc_east,phi_bc_south,phi_bc_north,phi_bc_low,phi_bc_high,\
1244 phi_bc_west_type,phi_bc_east_type,phi_bc_south_type,phi_bc_north_type,phi_bc_low_type,phi_bc_high_type,variable):
1245 if gpu:
1246 import cupy as xp
1247 else:
1248 import numpy as xp
1249
1250 phi3d_face_w=xp.empty((ni+1,nj,nk))
1251 phi3d_face_s=xp.empty((ni,nj+1,nk))
1252 phi3d_face_l=xp.empty((ni,nj,nk+1))
1253 phi3d_face_w[0:-1,:,:]=fx*phi3d+(1-fx)*xp.roll(phi3d,1,axis=0)
1254 phi3d_face_s[:,0:-1,:]=fy*phi3d+(1-fy)*xp.roll(phi3d,1,axis=1)
1255 phi3d_face_l[:,:,0:-1]=fz*phi3d+(1-fz)*xp.roll(phi3d,1,axis=2)
1256
1257
1258# west boundary
1259 phi3d_face_w[0,:,:]=phi_bc_west
1260 if phi_bc_west_type == 'n':
1261# neumann
1262 phi3d_face_w[0,:,:]=phi3d[0,:,:]
1263 if cyclic_x:
1264 phi3d_face_w[0,:,:]=0.5*(phi3d[0,:,:]+phi3d[-1,:,:])
1265
1266# east boundary
1267 phi3d_face_w[-1,:,:]=phi_bc_east
1268 if phi_bc_east_type == 'n':
1269# neumann
1270 phi3d_face_w[-1,:,:]=phi3d[-1,:,:]
1271 if cyclic_x:
1272 phi3d_face_w[-1,:,:]=0.5*(phi3d[0,:,:]+phi3d[-1,:,:])
1273
1274# south boundary
1275 phi3d_face_s[:,0,:]=phi_bc_south
1276 if phi_bc_south_type == 'n':
1277# neumann
1278 phi3d_face_s[:,0,:]=phi3d[:,0,:]
1279# d2phidy2=0
1280 if phi_bc_south_type == '2':
1281 phi3d_face_s[:,0,:]=1.5*phi3d[:,0,:]-0.5*phi3d[:,1,:]
1282
1283# north boundary
1284 phi3d_face_s[:,-1,:]=phi_bc_north
1285 if phi_bc_north_type == 'n':
1286# neumann
1287 phi3d_face_s[:,-1,:]=phi3d[:,-1,:]
1288 if phi_bc_north_type == '2':
1289# d2phidy2=0
1290 phi3d_face_s[:,-1,:]=1.5*phi3d[:,-1,:]-0.5*phi3d[:,-2,:]
1291
1292# low boundary
1293 phi3d_face_l[:,:,0]=phi_bc_low
1294# high boundary
1295 phi3d_face_l[:,:,-1]=phi_bc_high
1296 if phi_bc_low_type == 'n':
1297# neumann
1298# low boundary
1299 phi3d_face_l[:,:,0]= phi3d[:,:,0]
1300 if phi_bc_high_type == 'n':
1301# high boundary
1302 phi3d_face_l[:,:,-1]= phi3d[:,:,-1]
1303 if cyclic_z:
1304# low boundary
1305 phi3d_face_l[:,:,0]= 0.5*(phi3d[:,:,-1]+phi3d[:,:,0])
1306# high boundary
1307 phi3d_face_l[:,:,-1]= 0.5*(phi3d[:,:,-1]+phi3d[:,:,0])
1308
1309# this is needed only when blockage ios used
1310 if i_block_start !=0 or i_block_end !=0 or j_block_start !=0 or i_block_end !=0:
1311 phi3d_face_w,phi3d_face_s,phi3d_face_l = modify_face(phi3d_face_w,phi3d_face_s,phi3d_face_l,variable)
1312
1313 return phi3d_face_w,phi3d_face_s,phi3d_face_l
1314
1315def dphidx(phi_face_w,phi_face_s):
1316
1317 phi_w=phi_face_w[0:-1,:,:]*areawx[0:-1,:,:]
1318 phi_e=-phi_face_w[1:,:,:]*areawx[1:,:,:]
1319 phi_s=phi_face_s[:,0:-1,:]*areasx[:,0:-1,:]
1320 phi_n=-phi_face_s[:,1:,:]*areasx[:,1:,:]
1321 return (phi_w+phi_e+phi_s+phi_n)/vol
1322
1323def dphidy(phi_face_w,phi_face_s):
1324
1325 phi_w=phi_face_w[0:-1,:,:]*areawy[0:-1,:,:]
1326 phi_e=-phi_face_w[1:,:,:]*areawy[1:,:,:]
1327 phi_s=phi_face_s[:,0:-1,:]*areasy[:,0:-1,:]
1328 phi_n=-phi_face_s[:,1:,:]*areasy[:,1:,:]
1329 return (phi_w+phi_e+phi_s+phi_n)/vol
1330
1331def dphidz(phi_face_l):
1332
1333 phi_l=phi_face_l[:,:,0:-1]*areal[:,:,0:-1]
1334 phi_h=phi_face_l[:,:,1:]*areal[:,:,1:]
1335 return (phi_h-phi_l)/vol
1336
1337def coeff_m(convw,convs,convl,vis3d,u3d,v3d,w3d):
1338
1339 if itstep == 0 and iter == 0:
1340 print('muscle scheme')
1341
1342 visw=xp.zeros((ni+1,nj,nk))
1343 viss=xp.zeros((ni,nj+1,nk))
1344 visl=xp.zeros((ni,nj,nk+1))
1345 vis_turb=(vis3d-viscos)
1346
1347 visw[0:-1,:,:]=fx*vis_turb+(1-fx)*xp.roll(vis_turb,1,axis=0)+viscos
1348 viss[:,0:-1,:]=fy*vis_turb+(1-fy)*xp.roll(vis_turb,1,axis=1)+viscos
1349 visl[:,:,0:-1]=fz*vis_turb+(1-fz)*xp.roll(vis_turb,1,axis=2)+viscos
1350
1351 volw=xp.ones((ni+1,nj,nk))*1e-10
1352 vols=xp.ones((ni,nj+1,nk))*1e-10
1353 voll=xp.ones((ni,nj,nk+1))*1e-10
1354 volw[1:,:,:]=0.5*xp.roll(vol,-1,axis=0)+0.5*vol
1355 diffw=visw[0:-1,:,:]*areaw[0:-1,:,:]**2/volw[0:-1,:,:]
1356 vols[:,1:,:]=0.5*xp.roll(vol,-1,axis=1)+0.5*vol
1357 diffs=viss[:,0:-1,:]*areas[:,0:-1,:]**2/vols[:,0:-1,:]
1358 voll[:,:,1:]=0.5*xp.roll(vol,-1,axis=2)+0.5*vol
1359 diffl=visl[:,:,0:-1]*areal[:,:,0:-1]**2/voll[:,:,0:-1]
1360
1361 if cyclic_x:
1362 visw[0,:,:]=0.5*(vis_turb[0,:,:]+vis_turb[-1,:,:])+viscos
1363 diffw[0,:,:]=visw[0,:,:]*areaw[0,:,:]**2/(0.5*(vol[0,:,:]+vol[-1,:,:]))
1364 if cyclic_z:
1365 visl[:,:,0]=0.5*(vis_turb[:,:,0]+vis_turb[:,:,-1])+viscos
1366
1367# cep=0.5+sign(0.5,conve(i,j,k))
1368# cwp=0.5+sign(0.5,conve(i-1,j,k))
1369# cem=sign(0.5,conve(i,j,k))-0.5
1370
1371 cwp=0.5+0.5*xp.sign(convw[0:-1,:,:])
1372 cep=0.5+0.5*xp.sign(convw[1:,:,:])
1373 cwm=0.5*xp.sign(convw[0:-1,:,:])-0.5
1374 cem=0.5*xp.sign(convw[1:,:,:])-0.5
1375
1376 csp=0.5+0.5*xp.sign(convs[:,0:-1,:])
1377 cnp=0.5+0.5*xp.sign(convs[:,1:,:])
1378 csm=0.5*xp.sign(convs[:,0:-1,:])-0.5
1379 cnm=0.5*xp.sign(convs[:,1:,:])-0.5
1380
1381 clp=0.5+0.5*xp.sign(convl[:,:,0:-1])
1382 chp=0.5+0.5*xp.sign(convl[:,:,1:])
1383 clm=0.5*xp.sign(convl[:,:,0:-1])-0.5
1384 chm=0.5*xp.sign(convl[:,:,1:])-0.5
1385
1386# fix boundaries: no contribution
1387 if not cyclic_x:
1388 cem[-1,:,:]=0
1389 cwp[0,:,:]=0
1390 cnm[:,-1,:]=0
1391 csp[:,0,:]=0
1392 if not cyclic_z:
1393 chm[:,:,-1]=0
1394 clp[:,:,0]=0
1395
1396# first-order upwind in left-hand side
1397 aw3d=xp.maximum(convw[0:-1,:,:],0)+diffw
1398 ae3d=xp.maximum(-convw[1:,:,:],-0)+xp.roll(diffw,-1,axis=0)
1399 as3d=xp.maximum(convs[:,0:-1,:],0)+diffs
1400 an3d=xp.maximum(-convs[:,1:,:],0)+xp.roll(diffs,-1,axis=1)
1401 al3d=xp.maximum(convl[:,:,0:-1],0)+diffl
1402 ah3d=xp.maximum(-convl[:,:,1:],0)+xp.roll(diffl,-1,axis=2)
1403
1404 su3d =muscle_source(u3d,cep,cem,cwp,cwm,cnp,cnm,csp,csm,chp,chm,clp,clm)
1405 su3d_v=muscle_source(v3d,cep,cem,cwp,cwm,cnp,cnm,csp,csm,chp,chm,clp,clm)
1406 su3d_w=muscle_source(w3d,cep,cem,cwp,cwm,cnp,cnm,csp,csm,chp,chm,clp,clm)
1407
1408 su3d=(1-blend)*su3d
1409 su3d_v=(1-blend)*su3d_v
1410 su3d_w=(1-blend)*su3d_w
1411
1412 if blend > 0:
1413
1414# central differencing
1415 aw3d_c=diffw+(1-fx)*convw[0:-1,:,:]
1416 ae3d_c=xp.roll(diffw,-1,axis=0)-xp.roll(fx,-1,axis=0)*convw[1:,:,:]
1417
1418 as3d_c=diffs+(1-fy)*convs[:,0:-1,:]
1419 an3d_c=xp.roll(diffs,-1,axis=1)-xp.roll(fy,-1,axis=1)*convs[:,1:,:]
1420
1421 al3d_c=diffl+(1-fz)*convl[:,:,0:-1]
1422 ah3d_c=xp.roll(diffl,-1,axis=2)-xp.roll(fz,-1,axis=2)*convl[:,:,1:]
1423
1424 if not cyclic_x:
1425 aw3d_c[0,:,:]=0
1426 ae3d_c[-1,:,:]=0
1427 aw3d[0,:,:]=0
1428 ae3d[-1,:,:]=0
1429
1430 as3d_c[:,0,:]=0
1431 an3d_c[:,-1,:]=0
1432 as3d[:,0,:]=0
1433 an3d[:,-1,:]=0
1434
1435 if not cyclic_z:
1436 al3d_c[:,:,0]=0
1437 ah3d_c[:,:,-1]=0
1438 al3d[:,:,0]=0
1439 ah3d[:,:,-1]=0
1440
1441# blend of CDS and muscle. blend=1 means fully CDS, deferred
1442# su_u=c*(
1443# . (ae_c-ae(i,j,k))*(acr*(phi(i+1,j,k,n)-phi(i,j,k,n))
1444# . +(1.-acr)*(phio(i+1,j,k,n)-phio(i,j,k,n)))
1445# .+(aw_c-aw(i,j,k))*(acr*(phi(i-1,j,k,n)-phi(i,j,k,n))
1446# . +(1.-acr)*(phio(i-1,j,k,n)-phio(i,j,k,n)))
1447
1448 a=acrank_conv
1449
1450# u equation
1451 su3d=su3d+blend*\
1452 ((ae3d_c-ae3d)*(a*(xp.roll(u3d,-1,axis=0)-u3d)+(1-a)*(xp.roll(u3d_old,-1,axis=0)-u3d_old))
1453 +(aw3d_c-aw3d)*(a*(xp.roll(u3d, 1,axis=0)-u3d)+(1-a)*(xp.roll(u3d_old, 1,axis=0)-u3d_old))
1454 +(an3d_c-an3d)*(a*(xp.roll(u3d,-1,axis=1)-u3d)+(1-a)*(xp.roll(u3d_old,-1,axis=1)-u3d_old))
1455 +(as3d_c-as3d)*(a*(xp.roll(u3d, 1,axis=1)-u3d)+(1-a)*(xp.roll(u3d_old, 1,axis=1)-u3d_old))
1456 +(ah3d_c-ah3d)*(a*(xp.roll(u3d,-1,axis=2)-u3d)+(1-a)*(xp.roll(u3d_old,-1,axis=2)-u3d_old))
1457 +(al3d_c-al3d)*(a*(xp.roll(u3d, 1,axis=2)-u3d)+(1-a)*(xp.roll(u3d_old, 1,axis=2)-u3d_old)))
1458
1459# v equation
1460 su3d_v=su3d_v+blend* \
1461 ((ae3d_c-ae3d)*(a*(xp.roll(v3d,-1,axis=0)-v3d)+(1-a)*(xp.roll(v3d_old,-1,axis=0)-v3d_old))
1462 +(aw3d_c-aw3d)*(a*(xp.roll(v3d, 1,axis=0)-v3d)+(1-a)*(xp.roll(v3d_old, 1,axis=0)-v3d_old))
1463 +(an3d_c-an3d)*(a*(xp.roll(v3d,-1,axis=1)-v3d)+(1-a)*(xp.roll(v3d_old,-1,axis=1)-v3d_old))
1464 +(as3d_c-as3d)*(a*(xp.roll(v3d, 1,axis=1)-v3d)+(1-a)*(xp.roll(v3d_old, 1,axis=1)-v3d_old))
1465 +(ah3d_c-ah3d)*(a*(xp.roll(v3d,-1,axis=2)-v3d)+(1-a)*(xp.roll(v3d_old,-1,axis=2)-v3d_old))
1466 +(al3d_c-al3d)*(a*(xp.roll(v3d, 1,axis=2)-v3d)+(1-a)*(xp.roll(v3d_old, 1,axis=2)-v3d_old)))
1467
1468# w equation
1469 su3d_w=su3d_w+blend*\
1470 ((ae3d_c-ae3d)*(a*(xp.roll(w3d,-1,axis=0)-w3d)+(1-a)*(xp.roll(w3d_old,-1,axis=0)-w3d_old))
1471 +(aw3d_c-aw3d)*(a*(xp.roll(w3d, 1,axis=0)-w3d)+(1-a)*(xp.roll(w3d_old, 1,axis=0)-w3d_old))
1472 +(an3d_c-an3d)*(a*(xp.roll(w3d,-1,axis=1)-w3d)+(1-a)*(xp.roll(w3d_old,-1,axis=1)-w3d_old))
1473 +(as3d_c-as3d)*(a*(xp.roll(w3d, 1,axis=1)-w3d)+(1-a)*(xp.roll(w3d_old, 1,axis=1)-w3d_old))
1474 +(ah3d_c-ah3d)*(a*(xp.roll(w3d,-1,axis=2)-w3d)+(1-a)*(xp.roll(w3d_old,-1,axis=2)-w3d_old))
1475 +(al3d_c-al3d)*(a*(xp.roll(w3d, 1,axis=2)-w3d)+(1-a)*(xp.roll(w3d_old, 1,axis=2)-w3d_old)))
1476
1477# use first-order upstream upstream of x=0.2
1478 if embedded:
1479 i1 = (xp.abs(x_embed-xp2d[:,0])).argmin() # find index which closest
1480 if itstep == 0 and iter == 0:
1481 print('embedded: i1',i1)
1482# aw3d=diffw+(1-fx)*convw[0:-1,:,:]
1483# ae3d=xp.roll(diffw,-1,axis=0)-xp.roll(fx,-1,axis=0)*convw[1:,:,:]
1484
1485# as3d=diffs+(1-fy)*convs[:,0:-1,:]
1486# an3d=xp.roll(diffs,-1,axis=1)-xp.roll(fy,-1,axis=1)*convs[:,1:,:]
1487
1488# al3d=diffl+(1-fz)*convl[:,:,0:-1]
1489# ah3d=xp.roll(diffl,-1,axis=2)-xp.roll(fz,-1,axis=2)*convl[:,:,1:]
1490
1491 aw3d[i1:,:,:]=diffw[i1:,:,:]+(1-fx[i1:,:,:])*convw[i1:-1,:,:]
1492 ae3d[i1:,:,:]=xp.roll(diffw[i1:,:,:],-1,axis=0)-xp.roll(fx[i1:,:,:],-1,axis=0)*convw[i1+1:,:,:]
1493
1494 as3d[i1:,:,:]=diffs[i1:,:,:]+(1-fy[i1:,:,:])*convs[i1:,0:-1,:]
1495 an3d[i1:,:,:]=xp.roll(diffs[i1:,:,:],-1,axis=1)-xp.roll(fy[i1:,:,:],-1,axis=1)*convs[i1:,1:,:]
1496
1497 al3d[i1:,:,:]=diffl[i1:,:,:]+(1-fz[i1:,:,:])*convl[i1:,:,0:-1]
1498 ah3d[i1:,:,:]=xp.roll(diffl[i1:,:,:],-1,axis=2)-xp.roll(fz[i1:,:,:],-1,axis=2)*convl[i1:,:,1:]
1499
1500 su3d[i1:,:,:]=0
1501 su3d_v[i1:,:,:]=0
1502 su3d_w[i1:,:,:]=0
1503
1504# aw3d[0:i1,:,:]=xp.maximum(convw[0:i1,:,:],0)+diffw[0:i1,:,:]
1505# ae3d[1:i1,:,:]=xp.maximum(-convw[1:i1,:,:],-0)+xp.roll(diffw[1:i1,:,:],-1,axis=0)
1506# as3d[:i1,:,:]=xp.maximum(convs[:i1,0:-1,:],0)+diffs[:i1,:,:]
1507# an3d[:i1,:,:]=xp.maximum(-convs[:i1,1:,:],0)+xp.roll(diffs[:i1,:,:],-1,axis=1)
1508# al3d[:i1,:,:]=xp.maximum(convl[:i1,:,0:-1],0)+diffl[:i1,:,:]
1509# ah3d[:i1,:,:]=xp.maximum(-convl[:i1,:,1:],0)+xp.roll(diffl[:i1,:,:],-1,axis=2)
1510
1511 if itstep == 0 and iter == 0:
1512 print('emnedded scheme; i_embed=',i1)
1513
1514
1515 apo3d=vol/dt[itstep]
1516
1517 if not cyclic_x:
1518 aw3d[0,:,:]=0
1519 ae3d[-1,:,:]=0
1520 as3d[:,0,:]=0
1521 an3d[:,-1,:]=0
1522 if not cyclic_z:
1523 al3d[:,:,0]=0
1524 ah3d[:,:,-1]=0
1525
1526 return aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,su3d_v,su3d_w
1527
1528def minmo(a,b):
1529
1530# asign=sign(1.,a)
1531# rminmo=asign*max(0.,min(abs(a),b*asign))
1532
1533 asign=xp.sign(a)
1534
1535 return asign*xp.maximum(0,xp.minimum(abs(a),b*asign))
1536
1537def muscle_source(phi3d,cep,cem,cwp,cwm,cnp,cnm,csp,csm,chp,chm,clp,clm):
1538
1539 phip=xp.roll(phi3d,-1,axis=0)
1540 phim=xp.roll(phi3d,1,axis=0)
1541 phipp=xp.roll(phi3d,-2,axis=0)
1542 phimm=xp.roll(phi3d,2,axis=0)
1543
1544# su(i,j,k)=su(i,j,k)-0.5*
1545# & (conve(i,j,k)*cep*rminmo(phie-phip,phip-phiw)
1546# & -conve(i,j,k)*cem*rminmo(phie-phip,phiee-phie)
1547# & -conve(i-1,j,k)*cwp*rminmo(phip-phiw,phiw-phiww)
1548# & +conve(i-1,j,k)*cwm*rminmo(phip-phiw,phie-phip)
1549
1550
1551
1552 ss=-0.5*(convw[1:,:,:]*(cep*minmo(phip-phi3d,phi3d-phim)-cem*minmo(phip-phi3d,phipp-phip)) \
1553 -convw[0:-1,:,:]*(cwp*minmo(phi3d-phim,phim-phimm)-cwm*minmo(phi3d-phim,phip-phi3d)))
1554
1555 phip=xp.roll(phi3d,-1,axis=1)
1556 phim=xp.roll(phi3d,1,axis=1)
1557 phipp=xp.roll(phi3d,-2,axis=1)
1558 phimm=xp.roll(phi3d,2,axis=1)
1559
1560 ss=ss\
1561 -0.5*(convs[:,1:,:]*(cnp*minmo(phip-phi3d,phi3d-phim)-cnm*minmo(phip-phi3d,phipp-phip)) \
1562 -convs[:,0:-1,:]*(csp*minmo(phi3d-phim,phim-phimm)-csm*minmo(phi3d-phim,phip-phi3d)))
1563
1564 phip=xp.roll(phi3d,-1,axis=2)
1565 phim=xp.roll(phi3d,1,axis=2)
1566 phipp=xp.roll(phi3d,-2,axis=2)
1567 phimm=xp.roll(phi3d,2,axis=2)
1568
1569 ss=ss\
1570 -0.5*(convl[:,:,1:]*(chp*minmo(phip-phi3d,phi3d-phim)-chm*minmo(phip-phi3d,phipp-phip)) \
1571 -convl[:,:,0:-1]*(clp*minmo(phi3d-phim,phim-phimm)-clm*minmo(phi3d-phim,phip-phi3d)))
1572
1573 return ss
1574
1575def minmo(a,b):
1576
1577# asign=sign(1.,a)
1578# rminmo=asign*max(0.,min(abs(a),b*asign))
1579
1580 asign=xp.sign(a)
1581
1582 return asign*xp.maximum(0,xp.minimum(abs(a),b*asign))
1583
1584def coeff(convw,convs,convl,vis3d,prand_1,prand_2,f1_sst,scheme_local):
1585
1586 if prand_1 == prand_2:
1587 prand = xp.ones((ni,nj,nk))*prand_1
1588 else:
1589 prand = f1_sst*prand_1+(1-f1_sst)*prand_2
1590
1591 visw=xp.zeros((ni+1,nj,nk))
1592 viss=xp.zeros((ni,nj+1,nk))
1593 visl=xp.zeros((ni,nj,nk+1))
1594 if prand_1 > 0:
1595 vis_turb=(vis3d-viscos)/prand
1596 else:
1597 fk3d_local=xp.maximum(fk3d,fkmin_limit) #this limit is used only in the diffusion
1598 vis_turb=(vis3d-viscos)/xp.abs(prand)/fk3d_local**2
1599
1600 visw[0:-1,:,:]=fx*vis_turb+(1-fx)*xp.roll(vis_turb,1,axis=0)+viscos
1601 viss[:,0:-1,:]=fy*vis_turb+(1-fy)*xp.roll(vis_turb,1,axis=1)+viscos
1602 visl[:,:,0:-1]=fz*vis_turb+(1-fz)*xp.roll(vis_turb,1,axis=2)+viscos
1603
1604 volw=xp.ones((ni+1,nj,nk))*1e-10
1605 vols=xp.ones((ni,nj+1,nk))*1e-10
1606 voll=xp.ones((ni,nj,nk+1))*1e-10
1607 volw[1:,:,:]=0.5*xp.roll(vol,-1,axis=0)+0.5*vol
1608 diffw=visw[0:-1,:,:]*areaw[0:-1,:,:]**2/volw[0:-1,:,:]
1609 vols[:,1:,:]=0.5*xp.roll(vol,-1,axis=1)+0.5*vol
1610 diffs=viss[:,0:-1,:]*areas[:,0:-1,:]**2/vols[:,0:-1,:]
1611 voll[:,:,1:]=0.5*xp.roll(vol,-1,axis=2)+0.5*vol
1612 diffl=visl[:,:,0:-1]*areal[:,:,0:-1]**2/voll[:,:,0:-1]
1613
1614 if cyclic_x:
1615 visw[0,:,:]=0.5*(vis_turb[0,:,:]+vis_turb[-1,:,:])+viscos
1616 diffw[0,:,:]=visw[0,:,:]*areaw[0,:,:]**2/(0.5*(vol[0,:,:]+vol[-1,:,:]))
1617 if cyclic_z:
1618 visl[:,:,0]=0.5*(vis_turb[:,:,0]+vis_turb[:,:,-1])+viscos
1619 diffl[:,:,0]=visl[:,:,0]*areal[:,:,0]**2/(0.5*(vol[:,:,0]+vol[:,:,-1]))
1620
1621 if scheme_local == 'h':
1622 if itstep == 0 and iter == 0:
1623 print('hybrid scheme, prand=',prand_1)
1624
1625 aw3d=xp.maximum(convw[0:-1,:,:],diffw+(1-fx)*convw[0:-1,:,:])
1626 aw3d=xp.maximum(aw3d,0.)
1627
1628 ae3d=xp.maximum(-convw[1:,:,:],xp.roll(diffw,-1,axis=0)-xp.roll(fx,-1,axis=0)*convw[1:,:,:])
1629 ae3d=xp.maximum(ae3d,0.)
1630
1631 as3d=xp.maximum(convs[:,0:-1,:],diffs+(1-fy)*convs[:,0:-1,:])
1632 as3d=xp.maximum(as3d,0.)
1633
1634 an3d=xp.maximum(-convs[:,1:,:],xp.roll(diffs,-1,axis=1)-xp.roll(fy,-1,axis=1)*convs[:,1:,:])
1635 an3d=xp.maximum(an3d,0.)
1636
1637 al3d=xp.maximum(convl[:,:,0:-1],diffl+(1-fz)*convl[:,:,0:-1])
1638 al3d=xp.maximum(al3d,0.)
1639
1640 ah3d=xp.maximum(-convl[:,:,1:],xp.roll(diffl,-1,axis=2)-xp.roll(fz,-1,axis=2)*convl[:,:,1:])
1641 ah3d=xp.maximum(ah3d,0.)
1642
1643
1644 if scheme_local == 'u':
1645 if itstep == 0 and iter == 0:
1646 print('upwind scheme, prand=',prand_1)
1647
1648 aw3d=xp.maximum(convw[0:-1,:,:],0)+diffw
1649 ae3d=xp.maximum(-convw[1:,:,:],-0)+xp.roll(diffw,-1,axis=0)
1650 as3d=xp.maximum(convs[:,0:-1,:],0)+diffs
1651 an3d=xp.maximum(-convs[:,1:,:],0)+xp.roll(diffs,-1,axis=1)
1652 al3d=xp.maximum(convl[:,:,0:-1],0)+diffl
1653 ah3d=xp.maximum(-convl[:,:,1:],0)+xp.roll(diffl,-1,axis=2)
1654
1655 if scheme_local == 'c':
1656 if itstep == 0 and iter == 0:
1657 print('CDS scheme, prand=',prand_1)
1658 aw3d=diffw+(1-fx)*convw[0:-1,:,:]
1659 ae3d=xp.roll(diffw,-1,axis=0)-xp.roll(fx,-1,axis=0)*convw[1:,:,:]
1660
1661 as3d=diffs+(1-fy)*convs[:,0:-1,:]
1662 an3d=xp.roll(diffs,-1,axis=1)-xp.roll(fy,-1,axis=1)*convs[:,1:,:]
1663
1664 al3d=diffl+(1-fz)*convl[:,:,0:-1]
1665 ah3d=xp.roll(diffl,-1,axis=2)-xp.roll(fz,-1,axis=2)*convl[:,:,1:]
1666
1667
1668 apo3d=vol/dt[itstep]
1669
1670 if not cyclic_x:
1671 aw3d[0,:,:]=0
1672 ae3d[-1,:,:]=0
1673 as3d[:,0,:]=0
1674 an3d[:,-1,:]=0
1675 if not cyclic_z:
1676 al3d[:,:,0]=0
1677 ah3d[:,:,-1]=0
1678
1679 return aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d
1680
1681def bc(su3d,sp3d,phi_bc_west,phi_bc_east,phi_bc_south,phi_bc_north,phi_bc_low,phi_bc_high,\
1682 phi_bc_west_type,phi_bc_east_type,phi_bc_south_type,phi_bc_north_type,phi_bc_low_type,phi_bc_high_type):
1683
1684 su3d=xp.zeros((ni,nj,nk))
1685 sp3d=xp.zeros((ni,nj,nk))
1686
1687#south
1688 if phi_bc_south_type == 'd':
1689 sp3d[:,0,:]=sp3d[:,0,:]-viscos*as_bound
1690 su3d[:,0,:]=su3d[:,0,:]+viscos*as_bound*phi_bc_south
1691
1692#north
1693 if phi_bc_north_type == 'd':
1694 sp3d[:,-1,:]=sp3d[:,-1,:]-viscos*an_bound
1695 su3d[:,-1,:]=su3d[:,-1,:]+viscos*an_bound*phi_bc_north
1696
1697#west
1698 if phi_bc_west_type == 'd' and not cyclic_x:
1699 sp3d[0,:,:]=sp3d[0,:,:]-viscos*aw_bound
1700 su3d[0,:,:]=su3d[0,:,:]+viscos*aw_bound*phi_bc_west
1701#east
1702 if phi_bc_east_type == 'd' and not cyclic_x:
1703 sp3d[-1,:,:]=sp3d[-1,:,:]-viscos*ae_bound
1704 su3d[-1,:,:]=su3d[-1,:,:]+viscos*ae_bound*phi_bc_east
1705
1706#low
1707 if phi_bc_low_type == 'd' and not cyclic_z:
1708 sp3d[:,:,0]=sp3d[:,:,0]-viscos*al_bound
1709 su3d[:,:,0]=su3d[:,:,0]+viscos*al_bound*phi_bc_low
1710
1711#high
1712 if phi_bc_high_type == 'd' and not cyclic_z:
1713 sp3d[:,:,-1]=sp3d[:,:,-1]-viscos*ah_bound
1714 su3d[:,:,-1]=su3d[:,:,-1]+viscos*ah_bound*phi_bc_high
1715
1716 return su3d,sp3d
1717
1718def conv(u3d,v3d,w3d,p3d_face_w,p3d_face_s,p3d_face_l):
1719
1720 if itstep == 0 and iter == 0:
1721 print('conv called')
1722#compute convection
1723
1724 dtt=dt[itstep]*acrank
1725 u3d_star=u3d+dphidx(p3d_face_w,p3d_face_s)*dtt
1726 v3d_star=v3d+dphidy(p3d_face_w,p3d_face_s)*dtt
1727 w3d_star=w3d+dphidz(p3d_face_l)*dtt
1728
1729 u3d_face_w,u3d_face_s,u3d_face_l=compute_face_phi(u3d_star,u_bc_west,u_bc_east,u_bc_south,u_bc_north,u_bc_low,u_bc_high,\
1730 u_bc_west_type,u_bc_east_type,u_bc_south_type,u_bc_north_type,u_bc_low_type,u_bc_high_type,'u')
1731 v3d_face_w,v3d_face_s,v3d_face_l=compute_face_phi(v3d_star,v_bc_west,v_bc_east,v_bc_south,v_bc_north,v_bc_low,v_bc_high,\
1732 v_bc_west_type,v_bc_east_type,v_bc_south_type,v_bc_north_type,v_bc_low_type,v_bc_high_type,'v')
1733 w3d_face_w,w3d_face_s,w3d_face_l=compute_face_phi(w3d_star,w_bc_west,w_bc_east,w_bc_south,w_bc_north,w_bc_low,w_bc_high,\
1734 w_bc_west_type,w_bc_east_type,w_bc_south_type,w_bc_north_type,w_bc_low_type,w_bc_high_type,'w')
1735
1736 convw=-u3d_face_w*areawx-v3d_face_w*areawy
1737 convs=-u3d_face_s*areasx-v3d_face_s*areasy
1738 convl=w3d_face_l*areal
1739
1740 convw,convs,convl=modify_conv(convw,convs,convl)
1741
1742 return convw,convs,convl
1743
1744def solve_3d(phi3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,tol_conv,nmax,acrank_conv_local,solver_local,variable):
1745 if itstep == 0 and iter == 0:
1746 print('solve_3d called')
1747 print('nmax,acrank_conv_local',nmax,acrank_conv_local)
1748
1749 start_time_solver = time.time()
1750
1751 aw=aw3d.flatten()*acrank_conv_local
1752 ae=ae3d.flatten()*acrank_conv_local
1753 as1=as3d.flatten()*acrank_conv_local
1754 an=an3d.flatten()*acrank_conv_local
1755 al=al3d.flatten()*acrank_conv_local
1756 ah=ah3d.flatten()*acrank_conv_local
1757 ap=ap3d.flatten()
1758
1759 m=ni*nj*nk
1760
1761
1762
1763 if cyclic_x and cyclic_z:
1764 al_cyc=xp.zeros(m)
1765 al_cyc[0:-1:nk]= al[0:-1:nk]
1766 al[0:-1:nk]=0
1767 ah_cyc=xp.zeros(m)
1768 ah_cyc[nk-1::nk]=ah[nk-1::nk]
1769 ah[nk-1:-1:nk]=0
1770 ah[-1]=0
1771 A = sparse.diags([ap, -ah[:-1], -al[1:], -al_cyc, -ah_cyc[nk-1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
1772 [0, 1, -1, nk-1, -(nk-1), nk, -nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
1773 elif not cyclic_z and cyclic_x:
1774 A = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
1775 [0, 1, -1, nk,-nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
1776 elif cyclic_z and not cyclic_x:
1777 al_cyc=xp.zeros(m)
1778 al_cyc[0:-1:nk]= al[0:-1:nk]
1779 al[0:-1:nk]=0
1780 ah_cyc=xp.zeros(m)
1781 ah_cyc[nk-1::nk]=ah[nk-1::nk]
1782 ah[nk-1::nk]=0
1783 A= sparse.diags([ap,-ah[:-1],-al[1:],-al_cyc, -ah_cyc[nk-1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:]], \
1784 [0, 1, -1, nk-1, -(nk-1), nk, -nk, nk*nj, -nk*nj], format='csr')
1785 elif not cyclic_z and not cyclic_x:
1786 A = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:]], [0, 1, -1, nk, -nk, nk*nj, -nk*nj], format='csr')
1787
1788 su=su3d.flatten()
1789 phi=phi3d.flatten()
1790
1791 res_su=xp.linalg.norm(su)
1792 resid_orig=xp.linalg.norm(A*phi - su)
1793
1794 phi_org=phi
1795
1796# bicg (BIConjugate Gradient)
1797# bicgstab (BIConjugate Gradient STABilized)
1798# cg (Conjugate Gradient) - symmetric positive definite matrices only
1799# cgs (Conjugate Gradient Squared)
1800# gmres (Generalized Minimal RESidual)
1801# minres (MINimum RESidual)
1802# qmr (Quasi
1803 resid=xp.linalg.norm(A*phi - su,ord=norm_order)
1804 tol=tol_conv
1805 abs_tol=1e-10
1806 if tol_conv < 0:
1807# use absolute convergence criterium
1808 abs_tol =abs(tol_conv)*resid
1809 tol_conv=0
1810
1811 if solver_local == 'cgs':
1812 if itstep == 0 and iter == 0:
1813 print('solver in solve_3d: cgs')
1814 phi,info=linalg.cgs(A,su,x0=phi, atol=tol_conv, rtol=tol, maxiter=nmax) # good
1815 if solver_local == 'cg':
1816 if itstep == 0 and iter == 0:
1817 print('solver in solve_3d: cg')
1818 phi,info=linalg.cg(A,su,x0=phi, atol=abs_tol, rtol=tol, maxiter=nmax) # good
1819 if solver_local == 'gmres':
1820 if itstep == 0 and iter == 0:
1821 print('solver in solve_3d: gmres')
1822 phi,info=linalg.gmres(A,su,x0=phi, atol=abs_tol, rtol=tol, maxiter=nmax) # good
1823 if solver_local == 'qmr':
1824 if itstep == 0 and iter == 0:
1825 print('solver in solve_3d: qmr')
1826 phi,info=linalg.qmr(A,su,x0=phi, atol=abs_tol, rtol=tol, maxiter=nmax) # good
1827 if solver_local == 'lgmres':
1828 if itstep == 0 and iter == 0:
1829 print('solver in solve_3d: lgmres,tol,atol',tol,tol_conv)
1830 phi,info=linalg.lgmres(A,su,x0=phi, atol=abs_tol, rtol=tol, maxiter=nmax) # good
1831 if info > 0:
1832 print('warning in module solve_3d: convergence in sparse matrix solver not reached')
1833 print('solve_3d, into',info)
1834# compute residual without normalizing with |b|=|su3d|
1835 resid=xp.linalg.norm(A*phi - su,ord=norm_order)
1836
1837 delta_phi=xp.max(xp.abs(phi-phi_org))
1838
1839 phi3d=xp.reshape(phi,(ni,nj,nk))
1840 phi3d_org=xp.reshape(phi_org,(ni,nj,nk))
1841
1842
1843 print('variable=',variable)
1844 print(f"{'residual history in solve_3d. Variable '}{variable}{': initial residual: '} {resid_orig:.2e}{'final residual: ':>20}{resid:.2e}\
1845 {'delta_phi: ':>15}{delta_phi:.2e}")
1846
1847 print(f"{'time solve_3d: '}{time.time()-start_time_solver:.2e}")
1848
1849 return phi3d,resid
1850
1851def solve_pyamg(phi3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,tol_conv,acrank_conv_local):
1852
1853 if itstep == 0 and iter == 0:
1854 print('solve_pyamg called,tol_conv=',tol_conv,'acrank_conv_local=',acrank_conv_local)
1855 print('relation method_phi, amg_cycle_phi=',amg_relax_phi,amg_cycle_phi)
1856
1857 aw=aw3d.flatten()*acrank_conv_local
1858 ae=ae3d.flatten()*acrank_conv_local
1859 as1=as3d.flatten()*acrank_conv_local
1860 an=an3d.flatten()*acrank_conv_local
1861 al=al3d.flatten()*acrank_conv_local
1862 ah=ah3d.flatten()*acrank_conv_local
1863 ap=ap3d.flatten()
1864
1865
1866 m=ni*nj*nk
1867
1868
1869
1870 if cyclic_x and cyclic_z:
1871 if itstep == 0 and iter == 0:
1872 print('cyclic_x cyclic_z')
1873 al_cyc=xp.zeros(m)
1874 al_cyc[0:-1:nk]= al[0:-1:nk]
1875 al[0:-1:nk]=0
1876 ah_cyc=xp.zeros(m)
1877 ah_cyc[nk-1::nk]=ah[nk-1::nk]
1878 ah[nk-1:-1:nk]=0
1879 ah[-1]=0
1880 A = sparse.diags([ap, -ah[:-1], -al[1:], -al_cyc, -ah_cyc[nk-1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
1881 [0, 1, -1, nk-1, -(nk-1), nk, -nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
1882 elif not cyclic_z and cyclic_x:
1883 if itstep == 0 and iter == 0:
1884 print('cyclic_x and not cyclic_z')
1885 A = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
1886 [0, 1, -1, nk,-nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
1887 elif cyclic_z and not cyclic_x:
1888 al_cyc=xp.zeros(m)
1889 al_cyc[0:-1:nk]= al[0:-1:nk]
1890 al[0:-1:nk]=0
1891 ah_cyc=xp.zeros(m)
1892 ah_cyc[nk-1::nk]=ah[nk-1::nk]
1893# ah[nk-1:-1:nk]=0
1894 ah[nk-1::nk]=0
1895 A= sparse.diags([ap,-ah[:-1],-al[1:],-al_cyc, -ah_cyc[nk-1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:]], \
1896 [0, 1, -1, nk-1, -(nk-1), nk, -nk, nk*nj, -nk*nj], format='csr')
1897 elif not cyclic_z and not cyclic_x:
1898 if itstep == 0 and iter == 0:
1899 print('not cyclic_z and not cyclic_x')
1900 A = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:]], [0, 1, -1, nk, -nk, nk*nj, -nk*nj], format='csr')
1901
1902
1903
1904 App = pyamg.ruge_stuben_solver(A) # construct the multigrid hierarchy
1905# Ap = pyamg.classical.ruge_stuben_solver(Ap)
1906
1907 print('in solve_pyamg')
1908
1909 phi=phi3d.flatten()
1910 su=su3d.flatten()
1911 phi_org=phi
1912 res_amg = []
1913 if amg_relax_phi == 'default':
1914 phi = App.solve(su, tol=tol_conv, x0=phi, residuals=res_amg)
1915 elif amg_relax_phi == 'direct':
1916 phi = linalg.spsolve(A,su)
1917 res_amg=xp.zeros(1)
1918 else:
1919 phi = App.solve(su, tol=tol_conv, x0=phi,accel=amg_relax_phi,cycle=amg_cycle_phi, residuals=res_amg)
1920
1921# if amg_relax_phi != 'direct':
1922 print('Residual history in pyAMG', ["%0.4e" % i for i in res_amg])
1923
1924# index_phi=xp.argmax(xp.abs(phi-phi_org))
1925 delta_phi=xp.max(xp.abs(phi-phi_org))
1926
1927 print(f"{'Residual history in solve_pyAMG: delta_phi: ':>25}{delta_phi:.2e}")
1928
1929 resid=xp.linalg.norm(A*phi - su,ord=norm_order)
1930
1931 phi3d=xp.reshape(phi,(ni,nj,nk))
1932
1933 return phi3d,resid
1934
1935def solve_pyamgx(phi3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,tol_conv,acrank_conv_local,variable):
1936
1937 global A_x,b_x,x_x,rsc_x,cfg_x,A
1938
1939 start_time_solver = time.time()
1940
1941 if itstep == 0 and iter == 0:
1942 print('solve_pyamgx called,tol_conv=',tol_conv,'acrank_conv_local=',acrank_conv_local,'variable=',variable)
1943
1944 su=su3d.flatten()
1945 phi=phi3d.flatten()
1946
1947 if variable == 'p' or variable == 'u' or variable == 'k' or variable == 'eps' or variable == 'om' or \
1948 (variable == 'v' and coeff_v) or (variable == 'w' and coeff_w):
1949# the coefficient matrix is the same for u, v, w. It is assumed that sp3d=0 (e.g. there must be no symmetry b.c.)
1950
1951 aw=aw3d.flatten()*acrank_conv_local
1952 ae=ae3d.flatten()*acrank_conv_local
1953 as1=as3d.flatten()*acrank_conv_local
1954 an=an3d.flatten()*acrank_conv_local
1955 al=al3d.flatten()*acrank_conv_local
1956 ah=ah3d.flatten()*acrank_conv_local
1957 ap=ap3d.flatten()
1958
1959 m=ni*nj*nk
1960
1961 if cyclic_x and cyclic_z:
1962 print('cyclic_x cyclic_z')
1963 al_cyc=xp.zeros(m)
1964 al_cyc[0:-1:nk]= al[0:-1:nk]
1965 al[0:-1:nk]=0
1966 ah_cyc=xp.zeros(m)
1967 ah_cyc[nk-1::nk]=ah[nk-1::nk]
1968 ah[nk-1:-1:nk]=0
1969 ah[-1]=0
1970 A = sparse.diags([ap, -ah[:-1], -al[1:], -al_cyc, -ah_cyc[nk-1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
1971 [0, 1, -1, nk-1, -(nk-1), nk, -nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
1972 elif not cyclic_z and cyclic_x:
1973 print('cyclic_x and not cyclic_z')
1974 A = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
1975 [0, 1, -1, nk,-nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
1976 elif cyclic_z and not cyclic_x:
1977 al_cyc=xp.zeros(m)
1978 al_cyc[0:-1:nk]= al[0:-1:nk]
1979 al[0:-1:nk]=0
1980 ah_cyc=xp.zeros(m)
1981 ah_cyc[nk-1::nk]=ah[nk-1::nk]
1982# ah[nk-1:-1:nk]=0
1983 ah[nk-1::nk]=0
1984 A= sparse.diags([ap,-ah[:-1],-al[1:],-al_cyc, -ah_cyc[nk-1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:]], \
1985 [0, 1, -1, nk-1, -(nk-1), nk, -nk, nk*nj, -nk*nj], format='csr')
1986 elif not cyclic_z and not cyclic_x:
1987 print('not cyclic_z and not cyclic_x')
1988 A = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:]], [0, 1, -1, nk, -nk, nk*nj, -nk*nj], format='csr')
1989
1990# Create matrices and vectors:
1991 A_x = pyamgx.Matrix().create(rsc)
1992 b_x = pyamgx.Vector().create(rsc)
1993 x_x = pyamgx.Vector().create(rsc)
1994
1995 A_x.upload_CSR(A)
1996 b_x.upload(su) #upload pyamg and rhs
1997 x_x.upload(phi)
1998 solverx.setup(A_x)
1999 print('solvers setup')
2000
2001
2002 b_x.upload(su) #upload pyamg and rhs
2003 x_x.upload(phi)
2004
2005# solve system:
2006 solverx.solve(b_x, x_x)
2007
2008# Download solution
2009 if xp is numpy:
2010 x_x.download(phi)
2011 else:
2012 x_x.download_raw(phi.data)
2013
2014 if (variable == 'p' or variable == 'w' or variable == 'k' or variable == 'eps' or variable == 'om')\
2015 or (variable == 'u' and coeff_v):
2016 A_x.destroy()
2017 b_x.destroy()
2018 x_x.destroy()
2019
2020 phi3d=xp.reshape(phi,(ni,nj,nk))
2021
2022 resid=xp.linalg.norm(A*phi - su,ord=norm_order)
2023
2024 phi3d=xp.reshape(phi,(ni,nj,nk))
2025
2026 print(f"{'resid in solve_pyamgx: '}{resid:.2e}")
2027
2028 print(f"{'time solve_pyamgx: '}{time.time()-start_time_solver:.2e}")
2029
2030 return phi3d,resid
2031
2032def solve_p(phi3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,tol_conv):
2033
2034 global Ap_solve_p
2035
2036 start_time_solver = time.time()
2037
2038 if itstep == 0 and iter == 0:
2039 print('solve_p called')
2040 print('relation method=',amg_relax)
2041
2042 if iter == 0 and itstep == 0:
2043 print('A and M computed,tol_conv=',tol_conv)
2044# aw=xp.matrix.flatten(aw3d)
2045# ae=xp.matrix.flatten(ae3d)
2046# as1=xp.matrix.flatten(as3d)
2047# an=xp.matrix.flatten(an3d)
2048# al=xp.matrix.flatten(al3d)
2049# ah=xp.matrix.flatten(ah3d)
2050# ap=xp.matrix.flatten(ap3d)
2051
2052 aw=aw3d.flatten()
2053 ae=ae3d.flatten()
2054 as1=as3d.flatten()
2055 an=an3d.flatten()
2056 al=al3d.flatten()
2057 ah=ah3d.flatten()
2058 ap=ap3d.flatten()
2059
2060 m=ni*nj*nk
2061
2062 if cyclic_x and cyclic_z:
2063 print('cyclic_x cyclic_z')
2064 al_cyc=xp.zeros(m)
2065 al_cyc[0:-1:nk]= al[0:-1:nk]
2066 al[0:-1:nk]=0
2067 ah_cyc=xp.zeros(m)
2068 ah_cyc[nk-1::nk]=ah[nk-1::nk]
2069 ah[nk-1:-1:nk]=0
2070 ah[-1]=0
2071 Ap = sparse.diags([ap, -ah[:-1], -al[1:], -al_cyc, -ah_cyc[nk-1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
2072 [0, 1, -1, nk-1, -(nk-1), nk, -nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
2073 elif not cyclic_z and cyclic_x:
2074 print('cyclic_x and not cyclic_z')
2075 Ap = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
2076 [0, 1, -1, nk,-nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
2077 elif cyclic_z and not cyclic_x:
2078 print('cyclic_z and not cyclic_x')
2079 al_cyc=xp.zeros(m)
2080 al_cyc[0:-1:nk]= al[0:-1:nk]
2081 al[0:-1:nk]=0
2082 ah_cyc=xp.zeros(m)
2083 ah_cyc[nk-1::nk]=ah[nk-1::nk]
2084# ah[nk-1:-1:nk]=0
2085 ah[nk-1::nk]=0
2086 Ap= sparse.diags([ap,-ah[:-1],-al[1:],-al_cyc, -ah_cyc[nk-1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:]], \
2087 [0, 1, -1, nk-1, -(nk-1), nk, -nk, nk*nj, -nk*nj], format='csr')
2088 elif not cyclic_z and not cyclic_x:
2089 print('not cyclic_z and not cyclic_x')
2090 Ap = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:]], [0, 1, -1, nk, -nk, nk*nj, -nk*nj], format='csr')
2091
2092
2093
2094 Ap_solve_p = pyamg.ruge_stuben_solver(Ap) # construct the multigrid hierarchy
2095
2096 print('in solve_p')
2097
2098# phi=xp.matrix.flatten(phi3d)
2099# su=xp.matrix.flatten(su3d)
2100 phi=phi3d.flatten()
2101 su=su3d.flatten()
2102 res_amg = []
2103 if amg_relax == 'default':
2104 phi = Ap_solve_p.solve(su, tol=tol_conv, x0=phi, residuals=res_amg)
2105 else:
2106 phi = Ap_solve_p.solve(su, tol=tol_conv, x0=phi, residuals=res_amg, accel=amg_relax,cycle=amg_cycle)
2107
2108 print('Residual history in solve_p', ["%0.4e" % i for i in res_amg])
2109
2110 phi3d=xp.reshape(phi,(ni,nj,nk))
2111
2112 print(f"{'time solve_p: '}{time.time()-start_time_solver:.2e}")
2113
2114 return phi3d
2115
2116
2117def solve_px(phi3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,tol_conv):
2118
2119 global Ap,A_xp,b_xp,x_xp,rsc,cfg,solverx_p
2120
2121 if itstep == 0 and iter == 0:
2122 print('solve_px called')
2123
2124 phi=phi3d.flatten()
2125 su=su3d.flatten()
2126
2127 if iter == 0 and itstep == 0:
2128 print('A and M computed,tol_conv=',tol_conv)
2129
2130 aw=aw3d.flatten()
2131 ae=ae3d.flatten()
2132 as1=as3d.flatten()
2133 an=an3d.flatten()
2134 al=al3d.flatten()
2135 ah=ah3d.flatten()
2136 ap=ap3d.flatten()
2137
2138 m=ni*nj*nk
2139
2140 if cyclic_x and cyclic_z:
2141 print('cyclic_x cyclic_z')
2142 al_cyc=xp.zeros(m)
2143 al_cyc[0:-1:nk]= al[0:-1:nk]
2144 al[0:-1:nk]=0
2145 ah_cyc=xp.zeros(m)
2146 ah_cyc[nk-1::nk]=ah[nk-1::nk]
2147 ah[nk-1:-1:nk]=0
2148 ah[-1]=0
2149 Ap = sparse.diags([ap, -ah[:-1], -al[1:], -al_cyc, -ah_cyc[nk-1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
2150 [0, 1, -1, nk-1, -(nk-1), nk, -nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
2151 elif not cyclic_z and cyclic_x:
2152 print('cyclic_x and not cyclic_z')
2153 Ap = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
2154 [0, 1, -1, nk,-nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
2155 elif cyclic_z and not cyclic_x:
2156 al_cyc=xp.zeros(m)
2157 al_cyc[0:-1:nk]= al[0:-1:nk]
2158 al[0:-1:nk]=0
2159 ah_cyc=xp.zeros(m)
2160 ah_cyc[nk-1::nk]=ah[nk-1::nk]
2161# ah[nk-1:-1:nk]=0
2162 ah[nk-1::nk]=0
2163 Ap= sparse.diags([ap,-ah[:-1],-al[1:],-al_cyc, -ah_cyc[nk-1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:]], \
2164 [0, 1, -1, nk-1, -(nk-1), nk, -nk, nk*nj, -nk*nj], format='csr')
2165 elif not cyclic_z and not cyclic_x:
2166 print('not cyclic_z and not cyclic_x')
2167 Ap = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:]], [0, 1, -1, nk, -nk, nk*nj, -nk*nj], format='csr')
2168
2169# Create matrices and vectors:
2170 A_xp = pyamgx.Matrix().create(rsc)
2171 b_xp = pyamgx.Vector().create(rsc)
2172 x_xp = pyamgx.Vector().create(rsc)
2173
2174 solverx_p = pyamgx.Solver().create(rsc, cfg_p)
2175
2176 A_xp.upload_CSR(Ap) #upload pyamg poisson problem
2177 b_xp.upload(su) #upload pyamg and rhs
2178 x_xp.upload(phi)
2179 solverx_p.setup(A_xp)
2180
2181 b_xp.upload(su) #upload pyamg and rhs
2182 x_xp.upload(phi)
2183
2184# solve system:
2185 solverx_p.solve(b_xp, x_xp)
2186
2187# Download solution
2188 if xp is numpy:
2189 x_xp.download(phi)
2190 else:
2191 x_xp.download_raw(phi.data)
2192
2193 phi3d=xp.reshape(phi,(ni,nj,nk))
2194
2195
2196 resid=xp.linalg.norm(Ap*phi - su,ord=norm_order)
2197
2198 print(f"{'resid in solve_px: '}{resid:.2e}")
2199
2200 return phi3d
2201
2202def calcu(su3d,sp3d,dpdx_old,p3d_face_w,p3d_face_s):
2203 if itstep == 0 and iter == 0:
2204 print('calcu called')
2205# b.c., sources, coefficients
2206
2207# presssure gradient
2208 dpdx=acrank*dphidx(p3d_face_w,p3d_face_s)+(1-acrank)*dpdx_old
2209
2210 su3d=su3d-dpdx*vol
2211
2212# modify su & sp
2213 su3d,sp3d=modify_u(su3d,sp3d)
2214# unsteady term added in crank_nicol
2215
2216 return su3d,sp3d
2217
2218# find max index
2219# ind = xp.unravel_index(xp.argmax(yplus_south, axis=None), yplus_south.shape)
2220
2221
2222
2223def calcv(su3d,sp3d,dpdy_old,p3d_face_w,p3d_face_s):
2224 if itstep == 0 and iter == 0:
2225 print('calcv called')
2226# b.c., sources, coefficients
2227
2228# presssure gradient
2229 dpdy=acrank*dphidy(p3d_face_w,p3d_face_s)+(1-acrank)*dpdy_old
2230 su3d=su3d-dpdy*vol
2231
2232# modify su & sp
2233 su3d,sp3d=modify_v(su3d,sp3d)
2234
2235 return su3d,sp3d
2236
2237def compute_fk(k3d,eps3d,fk3d):
2238
2239 if itstep == 0 and iter == 0:
2240 print('compute_fk called')
2241
2242 if pans:
2243 L_t=k3d**1.5/eps3d
2244 psi=xp.maximum(1,L_t/(cdes*delta_max))
2245
2246 fk3d=xp.maximum(1.-(psi-1.)/(c_eps_2-c_eps_1),0)
2247# fk3d=xp.maximum(1.-(psi-1.)/(c_eps_2-c_eps_1),fkmin_limit)
2248
2249 if keps_des:
2250 rl=k3d**1.5/eps3d
2251 fk3d=xp.maximum(1.,rl/(0.67*delta_max))
2252 if jl0 < 0:
2253 jl=xp.abs(jl0)
2254 fk3d[:,0:jl,:]=1
2255
2256 fk3d=modify_fk(fk3d)
2257
2258 return fk3d
2259
2260def calck_kom(su3d,sp3d,k3d,om3d,vis3d,u3d_face_w,u3d_face_s,v3d_face_w,v3d_face_s,w3d_face_l):
2261# b.c., sources, coefficients
2262 if itstep == 0 and iter == 0:
2263 print('calck_kom called')
2264
2265# production term
2266 dudx=dphidx(u3d_face_w,u3d_face_s)
2267 dvdx=dphidx(v3d_face_w,v3d_face_s)
2268 dwdx=dphidx(w3d_face_w,w3d_face_s)
2269
2270 dudy=dphidy(u3d_face_w,u3d_face_s)
2271 dvdy=dphidy(v3d_face_w,v3d_face_s)
2272 dwdy=dphidy(w3d_face_w,w3d_face_s)
2273
2274 dudz=dphidz(u3d_face_l)
2275 dvdz=dphidz(v3d_face_l)
2276 dwdz=dphidz(w3d_face_l)
2277
2278 gen= (2.*(dudx**2+dvdy**2+dwdz**2)+(dudz+dwdx)**2+(dvdz+dwdy)**2+(dudy+dvdx)**2)
2279 vist=xp.maximum(vis3d-viscos,1e-10)
2280 su3d=su3d+vist*gen*vol
2281
2282
2283# dissipation term
2284 if sst or kom_des:
2285 sp3d=sp3d-fk3d*cmu*om3d*vol
2286 else:
2287 sp3d=sp3d-cmu*om3d*vol
2288
2289# modify su & sp
2290 su3d,sp3d,comm_term=modify_k(su3d,sp3d,gen)
2291
2292# unsteady term added in crank_nicol
2293
2294 return su3d,sp3d,gen,comm_term
2295
2296def calcom(su3d,sp3d,om3d,gen,comm_term):
2297 if itstep == 0 and iter == 0:
2298 print('calcom called')
2299
2300
2301#--------production term
2302 su3d=su3d+c_omega_1*gen*vol
2303
2304#--------dissipation term
2305 sp3d=sp3d-c_omega_2*om3d*vol
2306
2307# modify su & sp
2308 su3d,sp3d=modify_om(su3d,sp3d,comm_term)
2309
2310 return su3d,sp3d
2311
2312def calcom_sst(su3d,sp3d,om3d,gen,comm_term):
2313 if itstep == 0 and iter == 0:
2314 print('calcom_sst called')
2315
2316 k3d_face_w,k3d_face_s,k3d_face_l=compute_face_phi(k3d,k_bc_west,k_bc_east,k_bc_south,k_bc_north,k_bc_low,k_bc_high,\
2317 k_bc_west_type,k_bc_east_type,k_bc_south_type,k_bc_north_type,k_bc_low_type,k_bc_high_type,'k')
2318 om3d_face_w,om3d_face_s,om3d_face_l=compute_face_phi(om3d,om_bc_west,om_bc_east,om_bc_south,om_bc_north,om_bc_low,om_bc_high,\
2319 om_bc_west_type,om_bc_east_type,om_bc_south_type,om_bc_north_type,om_bc_low_type,om_bc_highlow_type,'om')
2320
2321#-------- cross term
2322 dkdx=dphidx(k3d_face_w,k3d_face_s)
2323 dkdy=dphidy(k3d_face_w,k3d_face_s)
2324 dkdz=dphidz(k3d_face_l)
2325
2326 domdx=dphidx(om3d_face_w,om3d_face_s)
2327 domdy=dphidy(om3d_face_w,om3d_face_s)
2328 domdz=dphidz(om3d_face_l)
2329
2330 if om_bc_south_type == 'd':
2331 domdy[:,0,:]=-2*om_bc_south/dist3d[:,0,:]
2332
2333#--- north wall
2334 if om_bc_north_type == 'd':
2335 domdy[:,-1,:]=2*om_bc_north/dist3d[:,-1,:]
2336
2337 crosv=dkdx*domdx+dkdy*domdy+dkdz*domdz
2338 cross_term=2*cr_sst*crosv/om3d
2339
2340# f1_sst
2341 term1=xp.maximum(cross_term,1.e-10)
2342 term1b=500*viscos/(om3d*dist3d**2)
2343 term2=xp.maximum(k3d**0.5/(cmu*om3d*dist3d),term1b)
2344 term3=4*cr_sst*k3d/(term1*dist3d**2)
2345 zeta=xp.minimum(term2,term3)
2346 f1_sst=xp.tanh(zeta**4)
2347
2348# f2_sst
2349 zeta=xp.maximum(2*k3d**0.5/(cmu*om3d*dist3d),term1b)
2350 f2_sst=xp.tanh(zeta**2)
2351
2352 cross_term=cross_term*(1-f1_sst)
2353
2354 sp3d=sp3d+xp.minimum(cross_term,0)/om3d*vol
2355 su3d=su3d+xp.maximum(cross_term,0)*vol
2356
2357#------- interpolate constants
2358 c_sst_1 = f1_sst*c_omega_1_sst_1+(1-f1_sst)*c_omega_1_sst_2
2359 c_sst_2 = f1_sst*c_omega_2_sst_1+(1-f1_sst)*c_omega_2_sst_2
2360
2361 su3d=su3d+c_sst_1*gen*vol
2362
2363#--------dissipation term
2364 sp3d=sp3d-c_sst_2*om3d*vol
2365
2366# modify su & sp
2367 su3d,sp3d=modify_om(su3d,sp3d,comm_term)
2368
2369 return su3d,sp3d,f1_sst,f2_sst
2370
2371def calck(su3d,sp3d,k3d,eps3d,vis3d,u3d_face_w,u3d_face_s,v3d_face_w,v3d_face_s,w3d_face_l):
2372# b.c., sources, coefficients
2373 if itstep == 0 and iter == 0:
2374 print('calck called')
2375
2376# production term
2377 dudx=dphidx(u3d_face_w,u3d_face_s)
2378 dvdx=dphidx(v3d_face_w,v3d_face_s)
2379 dwdx=dphidx(w3d_face_w,w3d_face_s)
2380
2381 dudy=dphidy(u3d_face_w,u3d_face_s)
2382 dvdy=dphidy(v3d_face_w,v3d_face_s)
2383 dwdy=dphidy(w3d_face_w,w3d_face_s)
2384
2385 dudz=dphidz(u3d_face_l)
2386 dvdz=dphidz(v3d_face_l)
2387 dwdz=dphidz(w3d_face_l)
2388
2389 gen= (2.*(dudx**2+dvdy**2+dwdz**2)+(dudz+dwdx)**2+(dvdz+dwdy)**2+(dudy+dvdx)**2)
2390 vist=xp.maximum(vis3d-viscos,1e-10)
2391 su3d=su3d+vist*gen*vol
2392
2393# dissipation term
2394 if keps_des:
2395 sp3d=sp3d-fk3d*eps3d/k3d*vol
2396 elif k_eq_les:
2397 delta_les=vol**0.3333
2398 sp3d=sp3d-c_eps*k3d**0.5/delta_les*vol
2399 else:
2400 sp3d=sp3d-eps3d/k3d*vol
2401
2402# modify su & sp
2403 su3d,sp3d,comm_term=modify_k(su3d,sp3d,gen)
2404
2405# unsteady term added in crank_nicol
2406
2407 return su3d,sp3d,gen,dudx,dudy
2408
2409def calceps(su3d,sp3d,k3d,eps3d,vis3d,gen):
2410 if itstep == 0 and iter == 0:
2411 print('calceps called')
2412
2413# b.c., sources, coefficients
2414 ueps=(eps3d*viscos)**0.25
2415 ystar=ueps*dist3d/viscos
2416 rt=k3d**2/eps3d/viscos
2417 fdampf2=((1.-xp.exp(-ystar/3.1))**2)*(1.-0.3*xp.exp(-(rt/6.5)**2))
2418 fmu3d=((1.-xp.exp(-ystar/14.))**2)*(1.+5./rt**0.75*xp.exp(-(rt/200.)**2))
2419 fmu3d=xp.minimum(fmu3d,1.)
2420
2421#--------production term
2422 su3d=su3d+c_eps_1*cmu*fmu3d*gen*k3d*vol
2423 if pans:
2424 c2u=c_eps_1+fk3d*(fdampf2*c_eps_2-c_eps_1)
2425 else:
2426 c2u=fdampf2*c_eps_2
2427
2428#--------dissipation term
2429 sp3d=sp3d-c2u*eps3d*vol/k3d
2430
2431# modify su & sp
2432 su3d,sp3d=modify_eps(su3d,sp3d)
2433
2434 return su3d,sp3d,fmu3d
2435
2436def calcw(su3d,sp3d,dpdz_old,p3d_face_l):
2437# b.c., sources, coefficients
2438 if itstep == 0 and iter == 0:
2439 print('calcw called')
2440
2441# presssure gradient
2442 dpdz=acrank*dphidz(p3d_face_l)+(1-acrank)*dpdz_old
2443 su3d=su3d-dpdz*vol
2444
2445# modify su & sp
2446 su3d,sp3d=modify_w(su3d,sp3d)
2447
2448# unsteady term added in crank_nicol
2449 return su3d,sp3d
2450
2451def calcp(convw,convs,convl):
2452 if itstep == 0 and iter == 0:
2453 print('calcp called')
2454# b.c., sources, coefficients
2455 volw=xp.ones((ni+1,nj,nk))*1e-10
2456 vols=xp.ones((ni,nj+1,nk))*1e-10
2457 voll=xp.ones((ni,nj,nk+1))*1e-10
2458 volw[1:,:,:]=0.5*xp.roll(vol,-1,axis=0)+0.5*vol
2459 aw3d=areaw[0:-1,:,:]**2/volw[0:-1,:,:]
2460 vols[:,1:,:]=0.5*xp.roll(vol,-1,axis=1)+0.5*vol
2461 as3d=areas[:,0:-1,:]**2/vols[:,0:-1,:]
2462 voll[:,:,1:]=0.5*xp.roll(vol,-1,axis=2)+0.5*vol
2463 al3d=areal[:,:,0:-1]**2/voll[:,:,0:-1]
2464
2465 ae3d=xp.roll(aw3d,-1,axis=0)
2466 an3d=xp.roll(as3d,-1,axis=1)
2467 ah3d=xp.roll(al3d,-1,axis=2)
2468
2469
2470 if cyclic_x:
2471 aw3d[0,:,:]=areaw[0,:,:]**2/(0.5*(vol[0,:,:]+vol[-1,:,:]))
2472 ae3d[-1,:,:]=aw3d[0,:,:]
2473 else:
2474 aw3d[0,:,:]=0
2475 ae3d[-1,:,:]=0
2476
2477 if cyclic_z:
2478 al3d[:,:,0]=areal[:,:,0]**2/(0.5*(vol[:,:,0]+vol[:,:,-1]))
2479 ah3d[:,:,-1]=al3d[:,:,0]
2480 else:
2481 al3d[:,:,0]=0
2482 ah3d[:,:,-1]=0
2483
2484 as3d[:,0,:]=0
2485 an3d[:,-1,:]=0
2486
2487 ap3d=aw3d+ae3d+as3d+an3d+al3d+ah3d
2488
2489# set p3d=0 in [0,0,0] to make it non-singular
2490 ap3d[0,0,0]=1e10
2491
2492 return aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d
2493
2494
2495def correct_conv(u3d,v3d,w3d,p3d,aw3d_p,as3d_p,al3d_p):
2496# correct convections
2497# create ghost cells at east & west boundaries with Neumann b.c.
2498 p3d_w=p3d
2499 p3d_s=p3d
2500 p3d_l=p3d
2501 dtt=dt[itstep]*acrank
2502#\\\\\\\\\\\\\ west face
2503# set zeros and put if before row 0
2504 p3d_w=xp.concatenate((xp.zeros((1,nj,nk)),p3d_w),axis=0)
2505 if cyclic_x:
2506 convw[1:-1,:,:]=convw[1:-1,:,:]+aw3d_p[1:,:,:]*(p3d_w[1:-1,:,:]-p3d[1:,:,:])*dtt
2507 convw[0,:,:]=convw[0,:,:]+aw3d_p[0,:,:]*(p3d[-1,:,:]-p3d[0,:,:])*dtt
2508 convw[-1,:,:]=convw[0,:,:]
2509 else:
2510 convw[0:-1,:,:]=convw[0:-1,:,:]+aw3d_p*(p3d_w[0:-1,:,:]-p3d)*dtt
2511
2512
2513#\\\\\\\\\\\\\ south face
2514# set zeros and put it before column 0
2515 p3d_s=xp.concatenate((xp.zeros((ni,1,nk)),p3d_s),axis=1)
2516 convs[:,0:-1,:]=convs[:,0:-1,:]+as3d_p*(p3d_s[:,0:-1,:]-p3d)*dtt
2517
2518#\\\\\\\\\\\\\ low face
2519# set zeros and put it before column 0
2520 p3d_l=xp.concatenate((xp.zeros((ni,nj,1)),p3d_l),axis=2)
2521 if cyclic_z:
2522 convl[:,:,1:-1]=convl[:,:,1:-1]+al3d_p[:,:,1:]*(p3d_l[:,:,1:-1]-p3d[:,:,1:])*dtt
2523 convl[:,:,-1]=convl[:,:,-1]+al3d_p[:,:,-1]*(p3d[:,:,-1]-p3d[:,:,0])*dtt
2524 convl[:,:,0]=convl[:,:,-1]
2525 else:
2526 convl[:,:,0:-1]=convl[:,:,0:-1]+al3d_p*(p3d_l[:,:,0:-1]-p3d)*dtt
2527# boundary
2528 convl[:,:,0]=w3d_face_l[:,:,0]*areal[:,:,0]
2529
2530# continuity error
2531 su3d=convw[0:-1,:,:]-convw[1:,:,:]\
2532 +convs[:,0:-1,:]-convs[:,1:,:]\
2533 +convl[:,:,0:-1]-convl[:,:,1:]
2534
2535 return convw,convs,convl,p3d,u3d,v3d,w3d,su3d
2536
2537
2538def update(u3d,v3d,w3d,k3d,eps3d,om3d,p3d_face_w,p3d_face_s,p3d_face_l):
2539 u3d_old=u3d
2540 v3d_old=v3d
2541 w3d_old=w3d
2542 k3d_old=k3d
2543 eps3d_old=eps3d
2544 om3d_old=om3d
2545 dpdx_old=dphidx(p3d_face_w,p3d_face_s)
2546 dpdy_old=dphidy(p3d_face_w,p3d_face_s)
2547 dpdz_old=dphidz(p3d_face_l)
2548
2549 return u3d_old,v3d_old,w3d_old,k3d_old,eps3d_old,om3d_old,dpdx_old,dpdy_old,dpdz_old
2550
2551def time_stats(u3d_mean,v3d_mean,w3d_mean,p3d_mean,k3d_mean,eps3d_mean,om3d_mean,uu3d_stress,vv3d_stress,ww3d_stress,uv3d_stress,uw3d_stress,vw3d_stress,\
2552 fk3d_mean,vis3d_mean,gen_mean):
2553
2554 global itstep_stats_counter
2555
2556
2557 itstep_stats_counter=itstep_stats_counter+1
2558
2559 print('time_stats called: itstep_stats_counter=',itstep_stats_counter)
2560 u3d_mean=u3d_mean+u3d
2561 v3d_mean=v3d_mean+v3d
2562 w3d_mean=w3d_mean+w3d
2563 p3d_mean=p3d_mean+p3d
2564 k3d_mean=k3d_mean+k3d
2565 fk3d_mean=fk3d_mean+fk3d
2566 om3d_mean=om3d_mean+om3d
2567 eps3d_mean=eps3d_mean+eps3d
2568 vis3d_mean=vis3d_mean+vis3d
2569 uu3d_stress=uu3d_stress+u3d**2
2570 vv3d_stress=vv3d_stress+v3d**2
2571 ww3d_stress=ww3d_stress+w3d**2
2572 uv3d_stress=uv3d_stress+u3d*v3d
2573 uw3d_stress=uw3d_stress+u3d*w3d
2574 vw3d_stress=vw3d_stress+v3d*w3d
2575 gen_mean=gen_mean+gen
2576
2577 return u3d_mean,v3d_mean,w3d_mean,p3d_mean,k3d_mean,eps3d_mean,om3d_mean,uu3d_stress,vv3d_stress,ww3d_stress,uv3d_stress,uw3d_stress,vw3d_stress,\
2578 fk3d_mean,vis3d_mean,gen_mean
2579
2580def crank_nicol(phi3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv_local):
2581 ap3d=aw3d+ae3d+as3d+an3d+al3d+ah3d
2582 su3d=su3d+(apo3d-(1-acrank_conv_local)*ap3d)*phi3d_old
2583 ap3d=apo3d+acrank_conv_local*ap3d-sp3d
2584 su3d=su3d+(1-acrank_conv_local)*\
2585 (ae3d*xp.roll(phi3d_old,-1,axis=0)+aw3d*xp.roll(phi3d_old,1,axis=0) \
2586 +an3d*xp.roll(phi3d_old,-1,axis=1)+as3d*xp.roll(phi3d_old,1,axis=1) \
2587 +ah3d*xp.roll(phi3d_old,-1,axis=2)+al3d*xp.roll(phi3d_old,1,axis=2))
2588 return ap3d,su3d
2589
2590def vist_kom(vis3d,k3d,om3d):
2591 if itstep == 0 and iter == 0:
2592 print('vist_kom called')
2593
2594 visold= vis3d
2595 vis3d= k3d/om3d+viscos
2596
2597# modify viscosity
2598 vis3d=modify_vis(vis3d)
2599
2600# under-relax viscosity
2601 vis3d= urfvis*vis3d+(1.-urfvis)*visold
2602
2603 return vis3d
2604
2605def vist_k_eq(vis3d,k3d):
2606 if itstep == 0 and iter == 0:
2607 print('vist_k_eq called')
2608
2609 visold= vis3d
2610 delta_les=vol**0.3333
2611 vis3d= cmu*k3d**0.5*delta_les+viscos
2612
2613# modify viscosity
2614 vis3d=modify_vis(vis3d)
2615
2616# under-relax viscosity
2617 vis3d= urfvis*vis3d+(1.-urfvis)*visold
2618
2619 return vis3d
2620
2621def vist_sst(vis3d,k3d,om3d,f2_sst,gen):
2622 if itstep == 0 and iter == 0:
2623 print('vist_sst called')
2624
2625 visold= vis3d
2626 a1=cmu**0.5
2627 denom=xp.maximum(a1*om3d,gen**0.5*f2_sst)
2628 vis3d= a1*k3d/denom+viscos
2629
2630# modify viscosity
2631 vis3d=modify_vis(vis3d)
2632
2633
2634# under-relax viscosity
2635 vis3d= urfvis*vis3d+(1.-urfvis)*visold
2636
2637
2638 return vis3d
2639
2640def vist_pans(vis3d,k3d,eps3d,fmu3d):
2641 if itstep == 0 and iter == 0:
2642 print('vist_pans called')
2643
2644 vis3d= cmu*fmu3d*k3d**2/eps3d+viscos
2645# under-relax viscosity
2646
2647# modify viscosity
2648 vis3d=modify_vis(vis3d)
2649
2650 visold= vis3d
2651 vis3d= urfvis*vis3d+(1.-urfvis)*visold
2652
2653
2654
2655 return vis3d
2656
2657def vist_smag(u3d_face_w,u3d_face_s,u3d_face_l,v3d_face_w,v3d_face_s,v3d_face_l,w3d_face_w,w3d_face_s,w3d_face_l,vis3d):
2658 if itstep == 0 and iter == 0:
2659 print('vist_smag called')
2660 dudx=dphidx(u3d_face_w,u3d_face_s)
2661 dvdx=dphidx(v3d_face_w,v3d_face_s)
2662 dwdx=dphidx(w3d_face_w,w3d_face_s)
2663
2664 dudy=dphidy(u3d_face_w,u3d_face_s)
2665 dvdy=dphidy(v3d_face_w,v3d_face_s)
2666 dwdy=dphidy(w3d_face_w,w3d_face_s)
2667
2668 dudz=dphidz(u3d_face_l)
2669 dvdz=dphidz(v3d_face_l)
2670 dwdz=dphidz(w3d_face_l)
2671
2672 gen= (2.*(dudx**2+dvdy**2+dwdz**2)+(dudz+dwdx)**2+(dvdz+dwdy)**2+(dudy+dvdx)**2)
2673
2674# RANS lengthscale
2675 rl_rans=0.41*dist3d
2676 rl_les=cmu*vol**0.3333333
2677 l_dist=0.15*dist3d
2678 l_max=0.15*delta_max
2679 dy=xp.diff(y2d[1:,:],axis=1)
2680# make it 3d
2681 dy=xp.repeat(dy[:,:,None],repeats=nk,axis=2)
2682
2683 l_temp=xp.maximum(l_dist,l_max)
2684 l_temp=xp.maximum(l_temp,dy)
2685 l_les=xp.minimum(l_temp,delta_max)
2686
2687# rl=xp.minimum(rl_rans,cmu*l_les)
2688
2689# IDDES
2690 alpha=0.25-dist3d/delta_max
2691 f_b= xp.minimum(2.*xp.exp(-9*alpha**2),1.)
2692 rl=f_b*rl_rans+(1-f_b)*cmu*l_les
2693
2694 visold= vis3d
2695 vis3d= rl**2*gen**0.5+viscos
2696
2697# modify viscosity
2698 vis3d=modify_vis(vis3d)
2699
2700# under-relax viscosity
2701 vis3d= urfvis*vis3d+(1.-urfvis)*visold
2702 return vis3d
2703
2705 scalar_names = ['pressure']
2706 scalar_variables = [p3d]
2707 if keps or pans or keps_des:
2708 scalar_names.append('turb_kin')
2709 scalar_names.append('epsilon')
2710 scalar_variables.append(k3d)
2711 scalar_variables.append(eps3d)
2712 if k_eq_les:
2713 scalar_names.append('turb_kin')
2714 scalar_variables.append(k3d)
2715 if sst or kom or kom_des:
2716 scalar_names.append('turb_kin')
2717 scalar_names.append('omega')
2718 scalar_variables.append(k3d)
2719 scalar_variables.append(om3d)
2720
2721 if save_vtk_movie:
2722 file_name = '%s.%d.vtk' % (vtk_file_name, itstep)
2723 else:
2724 file_name = '%s.vtk' % (vtk_file_name)
2725
2726 f = open(file_name,'w')
2727 f.write('# vtk DataFile Version 3.0\npyCALC-LES Data\nASCII\nDATASET STRUCTURED_GRID\n')
2728 f.write('DIMENSIONS %d %d %d\nPOINTS %d double\n' % (nk+1,nj+1,ni+1,(ni+1)*(nj+1)*(nk+1)))
2729 for i in range(ni+1):
2730 for j in range(nj+1):
2731 for k in range(nk+1):
2732 f.write('%.5f %.5f %.5f\n' % (x2d[i,j],y2d[i,j],dz*k))
2733 f.write('\nCELL_DATA %d\n' % (ni*nj*nk))
2734
2735 f.write('\nVECTORS velocity double\n')
2736 for i in range(ni):
2737 for j in range(nj):
2738 for k in range(nk):
2739 f.write('%.12e %.12e %.12e\n' % (u3d[i,j,k],v3d[i,j,k],w3d[i,j,k]))
2740
2741 for v in range(len(scalar_names)):
2742 var_name = scalar_names[v]
2743 var = scalar_variables[v]
2744 f.write('\nSCALARS %s double 1\nLOOKUP_TABLE default\n' % (var_name))
2745 for i in range(ni):
2746 for j in range(nj):
2747 for k in range(nk):
2748 f.write('%.10e\n' % (var[i,j,k]))
2749 f.close()
2750
2751 print('Flow state save into VTK format to file %s\n' % (file_name))
2752
2753def save_time_aver_data(u3d_mean,v3d_mean,w3d_mean,p3d_mean,eps3d_mean,om3d_mean,fk3d_mean,vis3d_mean,k3d_mean,uu3d_stress, \
2754 vv3d_stress,ww3d_stress,uv3d_stress,uw3d_stress,vw3d_stress,gen_mean):
2755
2756 print('save_time_aver_data called')
2757# save time-averaged data to disk
2758# if equi-distant mesh in z direction
2759 if save_average_z:
2760 umm=xp.mean(u3d_mean,axis=2)
2761 xp.save('u_averaged', xp.mean(u3d_mean,axis=2))
2762 xp.save('v_averaged', xp.mean(v3d_mean,axis=2))
2763 xp.save('w_averaged', xp.mean(w3d_mean,axis=2))
2764 xp.save('p_averaged', xp.mean(p3d_mean,axis=2))
2765 xp.save('k_averaged', xp.mean(k3d_mean,axis=2))
2766 xp.save('fk_averaged', xp.mean(fk3d_mean,axis=2))
2767 xp.save('k_averaged', xp.mean(k3d_mean,axis=2))
2768 xp.save('om_averaged', xp.mean(om3d_mean,axis=2))
2769 xp.save('vis_averaged', xp.mean(vis3d_mean,axis=2))
2770 xp.save('eps_averaged', xp.mean(eps3d_mean,axis=2))
2771 xp.save('gen_averaged', xp.mean(gen_mean,axis=2))
2772 xp.save('uu_stress', xp.mean(uu3d_stress,axis=2))
2773 xp.save('vv_stress', xp.mean(vv3d_stress,axis=2))
2774 xp.save('ww_stress', xp.mean(ww3d_stress,axis=2))
2775 xp.save('uv_stress', xp.mean(uv3d_stress,axis=2))
2776 xp.save('uw_stress', xp.mean(uw3d_stress,axis=2))
2777 xp.save('vw_stress', xp.mean(vw3d_stress,axis=2))
2778 xp.save('itstep',xp.array([xp.array(itstep_stats_counter),xp.array(nk),dz3d[0,0,0]]))
2779 print('itstep_stats_counter,nk,dz',itstep_stats_counter,nk,dz3d[0,0,0])
2780 print('data averaged in z')
2781 else:
2782 xp.save('u_averaged_3d', u3d_mean)
2783 xp.save('v_averaged_3d', v3d_mean)
2784 xp.save('w_averaged_3d', w3d_mean)
2785 xp.save('p_averaged_3d', p3d_mean)
2786 xp.save('vis_averaged_3d', vis3d_mean)
2787 if keps or kom_des or keps_des or kom or sst or k_eq_les:
2788 xp.save('k_averaged_3d', k3d_mean)
2789 if kom_des or kom:
2790 xp.save('om_averaged_3d', om3d_mean)
2791 if keps or keps_des or sst:
2792 xp.save('eps_averaged_3d', eps3d_mean)
2793 if kom_des or keps_des or sst:
2794 xp.save('fk_averaged_3d', fk3d_mean)
2795 xp.save('uu_stress_3d', uu3d_stress)
2796 xp.save('vv_stress_3d', vv3d_stress)
2797 xp.save('ww_stress_3d', ww3d_stress)
2798 xp.save('uv_stress_3d', uv3d_stress)
2799 xp.save('uw_stress_3d', uw3d_stress)
2800 xp.save('vw_stress_3d', vw3d_stress)
2801 xp.save('itstep',[itstep_stats_counter,nk,dz3d[0,0,0]])
2802 print('itstep_stats_counter,nk,dz',itstep_stats_counter,nk,dz3d[0,0,0])
2803 print('data not averaged in z')
2804 return
2805
2806 return
2807
2808def read_restart_data(u3d,v3d,w3d,p3d,k3d,eps3d,om3d):
2809
2810 print('read_restart_data called')
2811
2812 u3d=xp.load('u3d_saved.npy')
2813 v3d=xp.load('v3d_saved.npy')
2814 w3d=xp.load('w3d_saved.npy')
2815 p3d=xp.load('p3d_saved.npy')
2816 if keps or pans or keps_des or k_eq_les:
2817 k3d=xp.load('k3d_saved.npy')
2818 if keps or pans or keps_des:
2819 eps3d=xp.load('eps3d_saved.npy')
2820 if sst or kom or kom_des:
2821 k3d=xp.load('k3d_saved.npy')
2822 om3d=xp.load('om3d_saved.npy')
2823
2824
2825 return u3d,v3d,w3d,p3d,k3d,eps3d,om3d
2826
2827def save_data(u3d,v3d,w3d,p3d,k3d,eps3d,om3d):
2828
2829 print('save_data called')
2830
2831 xp.save('u3d_saved', u3d)
2832 xp.save('v3d_saved', v3d)
2833 xp.save('w3d_saved', w3d)
2834 xp.save('p3d_saved', p3d)
2835 if keps or pans or keps_des or k_eq_les:
2836 xp.save('k3d_saved', k3d)
2837 if keps or pans or keps_des:
2838 xp.save('eps3d_saved', eps3d)
2839 if sst or kom or kom_des:
2840 xp.save('k3d_saved', k3d)
2841 xp.save('om3d_saved', om3d)
2842
2843
2844 return
2845
2846def vist_wale(u3d_face_w,u3d_face_s,u3d_face_l,v3d_face_w,v3d_face_s,v3d_face_l,w3d_face_w,w3d_face_s,w3d_face_l,vis3d):
2847 if itstep == 0 and iter == 0:
2848 print('vist_wale called')
2849
2850 dudx=dphidx(u3d_face_w,u3d_face_s)
2851 dvdx=dphidx(v3d_face_w,v3d_face_s)
2852 dwdx=dphidx(w3d_face_w,w3d_face_s)
2853
2854 dudy=dphidy(u3d_face_w,u3d_face_s)
2855 dvdy=dphidy(v3d_face_w,v3d_face_s)
2856 dwdy=dphidy(w3d_face_w,w3d_face_s)
2857
2858 dudz=dphidz(u3d_face_l)
2859 dvdz=dphidz(v3d_face_l)
2860 dwdz=dphidz(w3d_face_l)
2861
2862 s11=dudx
2863 s12=0.5*(dudy+dvdx)
2864 s13=0.5*(dudz+dwdx)
2865
2866 s21=s12
2867 s22=dvdy
2868 s23=0.5*(dvdz+dwdy)
2869
2870 s31=s13
2871 s32=s23
2872 s33=dwdz
2873
2874 g11=dudx
2875 g12=dudy
2876 g13=dudz
2877
2878 g21=dvdx
2879 g22=dvdy
2880 g23=dvdz
2881
2882 g31=dwdx
2883 g32=dwdy
2884 g33=dwdz
2885
2886#square of g_ij = g_ik g_kj
2887 g11_2=g11*g11+g12*g21+g13*g31
2888 g12_2=g11*g12+g12*g22+g13*g32
2889 g13_2=g11*g13+g12*g23+g13*g33
2890
2891 g21_2=g21*g11+g22*g21+g23*g31
2892 g22_2=g21*g12+g22*g22+g23*g32
2893 g23_2=g21*g13+g22*g23+g23*g33
2894
2895 g31_2=g31*g11+g32*g21+g33*g31
2896 g32_2=g31*g12+g32*g22+g33*g32
2897 g33_2=g31*g13+g32*g23+g33*g33
2898
2899 gkk_2=(g11_2+g22_2+g33_2)/3.
2900
2901 sd11=g11_2-gkk_2
2902 sd12=0.5*(g12_2+g21_2)
2903 sd13=0.5*(g13_2+g31_2)
2904 sd21=sd12
2905 sd22=g22_2-gkk_2
2906 sd23=0.5*(g23_2+g32_2)
2907
2908 sd31=sd13
2909 sd32=sd23
2910 sd33=g33_2-gkk_2
2911
2912 sijsij=s11*s11+s12*s12+s13*s13+\
2913 s21*s21+s22*s22+s23*s23+\
2914 s31*s31+s32*s32+s33*s33
2915
2916 sdijsdij=sd11*sd11+sd12*sd12+sd13*sd13+\
2917 sd21*sd21+sd22*sd22+sd23*sd23+\
2918 sd31*sd31+sd32*sd32+sd33*sd33
2919
2920 cm=(10.6*0.1**2)**0.5
2921
2922 term1=sdijsdij**1.5/xp.maximum((sijsij**2.5+sdijsdij**1.25),1e-10)
2923
2924# RANS lengthscale
2925 rl_rans=0.41*dist3d
2926 rl_les=cmu*vol**0.3333333
2927 l_dist=0.15*dist3d
2928 l_max=0.15*delta_max
2929 dy=xp.diff(y2d[1:,:],axis=1)
2930# make it 3d
2931 dy=xp.repeat(dy[:,:,None],repeats=nk,axis=2)
2932
2933 l_temp=xp.maximum(l_dist,l_max)
2934 l_temp=xp.maximum(l_temp,dy)
2935 l_les=xp.minimum(l_temp,delta_max)
2936
2937 rl=xp.minimum(rl_rans,cm*l_les)
2938
2939
2940 visold= vis3d
2941 delta=vol**0.333333
2942# vis3d= (cm*delta)**2*term1+viscos
2943 vis3d= rl**2*term1+viscos
2944
2945# modify viscosity
2946 vis3d=modify_vis(vis3d)
2947
2948# under-relax viscosity
2949 vis3d= urfvis*vis3d+(1.-urfvis)*visold
2950 return vis3d
2951
2952
2953def main():
2954
2955 datax= xp.loadtxt("x2d.dat")
2956 x=datax[0:-1]
2957 ni=int(datax[-1])
2958 datay= xp.loadtxt("y2d.dat")
2959 y=datay[0:-1]
2960 nj=int(datay[-1])
2961
2962 x2d=xp.zeros((ni+1,nj+1))
2963 y2d=xp.zeros((ni+1,nj+1))
2964
2965 x2d=xp.reshape(x,(ni+1,nj+1))
2966 y2d=xp.reshape(y,(ni+1,nj+1))
2967
2968 dataz=xp.loadtxt('z.dat')
2969 if dataz[1] < 0:
2970# non-equi-distans grid in z
2971 nk=xp.abs(dataz[1])
2972 nk=int(nk)
2973 z=dataz[2:]
2974 zmax=dataz[-1]
2975 dz=xp.diff(z)
2976# dz3d
2977 dz3d=xp.ones((ni+1,nj+1,nk))
2978# make it 2D
2979 dz3d=xp.repeat(dz[None,:], repeats=nj+1, axis=0)
2980# make it 3D
2981 dz3d=xp.repeat(dz3d[None,:,:], repeats=ni+1, axis=0)
2982 else:
2983# equi-distans grid in z
2984 zmax=dataz[0]
2985 nk=dataz[1]
2986 nk=int(nk)
2987 dz=zmax/nk
2988 dz3d=xp.ones((ni+1,nj+1,nk))*dz
2989 z = xp.linspace(0, zmax, nk+1)
2990
2991# compute cell centers
2992 xp2d=0.25*(x2d[0:-1,0:-1]+x2d[0:-1,1:]+x2d[1:,0:-1]+x2d[1:,1:])
2993 yp2d=0.25*(y2d[0:-1,0:-1]+y2d[0:-1,1:]+y2d[1:,0:-1]+y2d[1:,1:])
2994 zp=0.5*(z[0:-1]+z[1:])
2995
2996
2997# initialize geometric arrays
2998
2999 vol=xp.zeros((ni,nj,nk))
3000 areas=xp.zeros((ni,nj+1,nk))
3001 areasx=xp.zeros((ni,nj+1,nk))
3002 areasy=xp.zeros((ni,nj+1,nk))
3003 areaw=xp.zeros((ni+1,nj,nk))
3004 areawx=xp.zeros((ni+1,nj,nk))
3005 areawy=xp.zeros((ni+1,nj,nk))
3006 areal=xp.zeros((ni,nj,nk+1))
3007 as_bound=xp.zeros((ni,nk))
3008 an_bound=xp.zeros((ni,nk))
3009 aw_bound=xp.zeros((nj,nk))
3010 ae_bound=xp.zeros((nj,nk))
3011 al_bound=xp.zeros((ni,nj))
3012 ah_bound=xp.zeros((ni,nj))
3013 fx=xp.zeros((ni,nj,nk))
3014 fy=xp.zeros((ni,nj,nk))
3015 fz=xp.zeros((ni,nj,nk))
3016
3017 for itstep in range(0,ntstep):
3018
3019
3020 for iter in range(0,maxit):
3021
3022 start_time_iter = time.time()
3023# coefficients for velocities
3024 start_time = time.time()
3025# conpute inlet fluc
3026 if iter == 0 and not cyclic_x:
3027 u_bc_west,v_bc_west,w_bc_west,k_bc_west,eps_bc_west,om_bc_west,u3d_face_w,convw = modify_inlet()
3028# boundary conditions for u3d
3029 su3d,sp3d=bc(su3d,sp3d,u_bc_west,u_bc_east,u_bc_south,u_bc_north,u_bc_low,u_bc_high, \
3030 u_bc_west_type,u_bc_east_type,u_bc_south_type,u_bc_north_type,u_bc_low_type,u_bc_high_type)
3031 if scheme == 'm':
3032 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d_local,su3d_v,su3d_w=coeff_m(convw,convs,convl,vis3d,u3d,v3d,w3d)
3033 su3d=su3d+su3d_local
3034 else:
3035 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d=coeff(convw,convs,convl,vis3d,1,1,f1_sst,scheme)
3036
3037# u3d
3038 su3d,sp3d=calcu(su3d,sp3d,dpdx_old,p3d_face_w,p3d_face_s)
3039 ap3d,su3d=crank_nicol(u3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv)
3040
3041# this is needed only when blockage ios used
3042 if i_block_start !=0 or i_block_end !=0 or j_block_start !=0 or i_block_end !=0:
3043 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_block('u')
3044
3045 if solver_vel == 'pyamg':
3046 u3d,residual_u=solve_pyamg(u3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_u,acrank_conv)
3047 elif solver_vel == 'pyamgx':
3048 u3d,residual_u=solve_pyamgx(u3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_u,acrank_conv,'u')
3049 else:
3050 u3d,residual_u=solve_3d(u3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_u,nsweep_vel,acrank_conv,solver_vel,'u')
3051 print(f"{'time u: '}{time.time()-start_time:.2e}")
3052 start_time = time.time()
3053# v3d
3054# boundary conditions for v3d
3055 su3d,sp3d=bc(su3d,sp3d,v_bc_west,v_bc_east,v_bc_south,v_bc_north,v_bc_low,v_bc_high, \
3056 v_bc_west_type,v_bc_east_type,v_bc_south_type,v_bc_north_type,v_bc_low_type,v_bc_high_type)
3057 if scheme == 'm':
3058 su3d=su3d+su3d_v
3059 su3d,sp3d=calcv(su3d,sp3d,dpdy_old,p3d_face_w,p3d_face_s)
3060 ap3d,su3d=crank_nicol(v3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv)
3061
3062# this is needed only when blockage ios used
3063 if i_block_start !=0 or i_block_end !=0 or j_block_start !=0 or i_block_end !=0:
3064 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_block('v')
3065
3066 if solver_vel == 'pyamg':
3067 v3d,residual_v=solve_pyamg(v3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_v,acrank_conv)
3068 elif solver_vel == 'pyamgx':
3069 v3d,residual_v=solve_pyamgx(v3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_v,acrank_conv,'v')
3070 else:
3071 v3d,residual_v=solve_3d(v3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_v,nsweep_vel,acrank_conv,solver_vel,'v')
3072 print(f"{'time v: '}{time.time()-start_time:.2e}")
3073
3074 start_time = time.time()
3075# w3d
3076# boundary conditions for w3d
3077 su3d,sp3d=bc(su3d,sp3d,w_bc_west,w_bc_east,w_bc_south,w_bc_north,w_bc_low,w_bc_high, \
3078 w_bc_west_type,w_bc_east_type,w_bc_south_type,w_bc_north_type,w_bc_low_type,w_bc_high_type)
3079 if scheme == 'm':
3080 su3d=su3d+su3d_w
3081 su3d,sp3d=calcw(su3d,sp3d,dpdz_old,p3d_face_l)
3082 ap3d,su3d=crank_nicol(w3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv)
3083
3084# this is needed only when blockage ios used
3085 if i_block_start !=0 or i_block_end !=0 or j_block_start !=0 or i_block_end !=0:
3086 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_block('w')
3087
3088 if solver_vel == 'pyamg':
3089 w3d,residual_w=solve_pyamg(w3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_w,acrank_conv)
3090 elif solver_vel == 'pyamgx':
3091 w3d,residual_w=solve_pyamgx(w3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_w,acrank_conv,'w')
3092 else:
3093 w3d,residual_w=solve_3d(w3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_w,nsweep_vel,acrank_conv,solver_vel,'w')
3094 print(f"{'time w: '}{time.time()-start_time:.2e}")
3095
3096 start_time = time.time()
3097# p3d
3098 convw,convs,convl=conv(u3d,v3d,w3d,p3d_face_w,p3d_face_s,p3d_face_l)
3099
3100 if not cyclic_x:
3101 convw,u_bc_east = modify_outlet(convw)
3102
3103
3104# RHS
3105# continuity error
3106 su3d=(convw[0:-1,:,:]-convw[1:,:,:]\
3107 +convs[:,0:-1,:]-convs[:,1:,:]\
3108 +convl[:,:,0:-1]-convl[:,:,1:])/acrank/dt[itstep]
3109
3110#
3111 if iter == 0 and itstep == 0:
3112 aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d=calcp(convw,convs,convl)
3113# this is needed only when blockage ios used
3114 if i_block_start !=0 or i_block_end !=0 or j_block_start !=0 or i_block_end !=0:
3115 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_block('p')
3116# coefficients are (maybe) modified in 'fix_block': re-compute ap3d
3117 ap3d=aw3d+ae3d+as3d+an3d+al3d+ah3d
3118
3119 aw3d_p=aw3d
3120 as3d_p=as3d
3121 al3d_p=al3d
3122
3123 print('xp.sum(su3d)',xp.sum(su3d))
3124
3125 if solver_p == 'pyamg':
3126 p3d=solve_p(p3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_p)
3127 elif solver_p == 'pyamgx':
3128# the A matrix is re-computed in solve_pyamgx every time
3129 aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d=calcp(convw,convs,convl)
3130 p3d,dummy=solve_pyamgx(p3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_p,1,'p')
3131
3132 elif solver_p == 'pyamgx_p':
3133# the A matrix is computed only once (requires more GPUY memory)
3134 p3d=solve_px(p3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_p)
3135 else:
3136 print('solver_p = ',solver_p)
3137 print('no such solver')
3138 sys.exit()
3139
3140# correct u, v, w, p
3141 convw,convs,convl,p3d,u3d,v3d,w3d,su3d= correct_conv(u3d,v3d,w3d,p3d,aw3d_p,as3d_p,al3d_p)
3142
3143 res_1d=su3d.flatten()
3144 residual_p=xp.linalg.norm(res_1d,ord=1)
3145
3146 u3d_face_w,u3d_face_s,u3d_face_l=compute_face_phi(u3d,u_bc_west,u_bc_east,u_bc_south,u_bc_north,u_bc_low,u_bc_high,\
3147 u_bc_west_type,u_bc_east_type,u_bc_south_type,u_bc_north_type,u_bc_low_type,u_bc_high_type,'u')
3148 v3d_face_w,v3d_face_s,v3d_face_l=compute_face_phi(v3d,v_bc_west,v_bc_east,v_bc_south,v_bc_north,v_bc_low,v_bc_high,\
3149 v_bc_west_type,v_bc_east_type,v_bc_south_type,v_bc_north_type,v_bc_low_type,v_bc_high_type,'v')
3150 w3d_face_w,w3d_face_s,w3d_face_l=compute_face_phi(w3d,w_bc_west,w_bc_east,w_bc_south,w_bc_north,w_bc_low,w_bc_high,\
3151 w_bc_west_type,w_bc_east_type,w_bc_south_type,w_bc_north_type,w_bc_low_type,w_bc_high_type,'w')
3152 p3d_face_w,p3d_face_s,p3d_face_l=compute_face_phi(p3d,p_bc_west,p_bc_east,p_bc_south,p_bc_north,p_bc_low,p_bc_high,\
3153 p_bc_west_type,p_bc_east_type,p_bc_south_type,p_bc_north_type,p_bc_low_type,p_bc_high_type,'p')
3154
3155
3156 print(f"{'time p: '}{time.time()-start_time:.2e}")
3157
3158 start_time = time.time()
3159
3160 if sst or kom or kom_des:
3161 if sst:
3162 vis3d=vist_sst(vis3d,k3d,om3d,f2_sst,gen)
3163 else:
3164 vis3d=vist_kom(vis3d,k3d,om3d)
3165# boundary conditions for k3d
3166 su3d,sp3d=bc(su3d,sp3d,k_bc_west,k_bc_east,k_bc_south,k_bc_north,k_bc_low,k_bc_high, \
3167 k_bc_west_type,k_bc_east_type,k_bc_south_type,k_bc_north_type,k_bc_low_type,k_bc_high_type)
3168# coefficients
3169 if sst:
3170 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d=coeff(convw,convs,convl,vis3d,prand_k_sst_1,prand_k_sst_2,f1_sst,scheme_turb)
3171 else:
3172 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d=coeff(convw,convs,convl,vis3d,prand_k,prand_k,f1_sst,scheme_turb)
3173# k
3174 su3d,sp3d,gen,comm_term=calck_kom(su3d,sp3d,k3d,om3d,vis3d,u3d_face_w,u3d_face_s,v3d_face_w,v3d_face_s,w3d_face_l)
3175
3176
3177 ap3d,su3d=crank_nicol(k3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv_kom)
3178
3179# this is needed only when blockage ios used
3180 if i_block_start !=0 or i_block_end !=0 or j_block_start !=0 or i_block_end !=0:
3181 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_block('k')
3182
3183 if solver_turb == 'pyamg':
3184 k3d,residual_k=solve_pyamg(k3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_k,acrank_conv_kom)
3185 elif solver_turb == 'pyamgx':
3186 k3d,residual_k=solve_pyamgx(k3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_k,acrank_conv_kom,'k')
3187 else:
3188 k3d,residual_k=solve_3d(k3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_k,nsweep_kom,acrank_conv_kom,solver_turb,'k')
3189 print(f"{'time k: '}{time.time()-start_time:.2e}")
3190
3191 k3d=xp.maximum(k3d,k_min)
3192
3193 start_time = time.time()
3194# omega
3195# boundary conditions for om3d
3196 su3d,sp3d=bc(su3d,sp3d,om_bc_west,om_bc_east,om_bc_south,om_bc_north,om_bc_low,om_bc_high, \
3197 om_bc_west_type,om_bc_east_type,om_bc_south_type,om_bc_north_type,om_bc_low_type,om_bc_high_type)
3198 if sst:
3199 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d=coeff(convw,convs,convl,vis3d,prand_omega_sst_1,prand_omega_sst_2,f1_sst,scheme_turb)
3200 else:
3201 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d=coeff(convw,convs,convl,vis3d,prand_omega,prand_omega,f1_sst,scheme_turb)
3202 if sst:
3203 su3d,sp3d,f1_sst,f2_sst= calcom_sst(su3d,sp3d,om3d,gen,comm_term)
3204 else:
3205 su3d,sp3d= calcom(su3d,sp3d,om3d,gen,comm_term)
3206 ap3d,su3d=crank_nicol(om3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv_kom)
3207
3208# this is needed only when blockage ios used
3209 if i_block_start !=0 or i_block_end !=0 or j_block_start !=0 or i_block_end !=0:
3210 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_block('om')
3211
3212 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_omega()
3213
3214 if solver_turb == 'pyamg':
3215 om3d,residual_om=solve_pyamg(om3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_om,acrank_conv_kom)
3216 elif solver_turb == 'pyamgx':
3217 om3d,residual_om=solve_pyamgx(om3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_om,acrank_conv_kom,'om')
3218 else:
3219 om3d,residual_om=solve_3d(om3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_om,nsweep_kom,acrank_conv_kom,solver_turb,'om')
3220 om3d=xp.maximum(om3d,om_min)
3221
3222 print(f"{'time omega: '}{time.time()-start_time:.2e}")
3223
3224 start_time = time.time()
3225
3226 if smag:
3227 vis3d=vist_smag(u3d_face_w,u3d_face_s,u3d_face_l,v3d_face_w,v3d_face_s,v3d_face_l,w3d_face_w,w3d_face_s,w3d_face_l,vis3d)
3228 if wale:
3229 vis3d=vist_wale(u3d_face_w,u3d_face_s,u3d_face_l,v3d_face_w,v3d_face_s,v3d_face_l,w3d_face_w,w3d_face_s,w3d_face_l,vis3d)
3230
3231 print(f"{'time Smag or Wale: '}{time.time()-start_time:.2e}")
3232 start_time = time.time()
3233 if pans or keps or keps_des or k_eq_les:
3234 if k_eq_les:
3235 vis3d=vist_k_eq(vis3d,k3d)
3236 else:
3237 vis3d=vist_pans(vis3d,k3d,eps3d,fmu3d)
3238# coefficients
3239 start_time = time.time()
3240 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d=coeff(convw,convs,convl,vis3d,prand_k,prand_k,f1_sst,scheme_turb)
3241# k
3242# boundary conditions for k3d
3243 su3d,sp3d=bc(su3d,sp3d,k_bc_west,k_bc_east,k_bc_south,k_bc_north,k_bc_low,k_bc_high, \
3244 k_bc_west_type,k_bc_east_type,k_bc_south_type,k_bc_north_type,k_bc_low_type,k_bc_high_type)
3245 su3d,sp3d,gen,dudx,dudy=calck(su3d,sp3d,k3d,eps3d,vis3d,u3d_face_w,u3d_face_s,v3d_face_w,v3d_face_s,w3d_face_l)
3246
3247 ap3d,su3d=crank_nicol(k3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv_keps)
3248
3249# this is needed only when blockage ios used
3250 if i_block_start !=0 or i_block_end !=0 or j_block_start !=0 or i_block_end !=0:
3251 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_block('k')
3252
3253 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_k()
3254
3255 if solver_turb == 'pyamg':
3256 k3d,residual_k=solve_pyamg(k3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_k,acrank_conv_keps)
3257 elif solver_turb == 'pyamgx':
3258 k3d,residual_k=solve_pyamgx(k3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_k,acrank_conv_keps,'k')
3259 else:
3260 k3d,residual_k=solve_3d(k3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_k,nsweep_keps,acrank_conv_keps,solver_turb,'k')
3261
3262
3263 k3d=xp.maximum(k3d,k_min)
3264 print(f"{'time k: '}{time.time()-start_time:.2e}")
3265 start_time = time.time()
3266
3267 if not k_eq_les:
3268 start_time = time.time()
3269# boundary conditions for eps3d
3270 su3d,sp3d=bc(su3d,sp3d,eps_bc_west,eps_bc_east,eps_bc_south,eps_bc_north,eps_bc_low,eps_bc_high, \
3271 eps_bc_west_type,eps_bc_east_type,eps_bc_south_type,eps_bc_north_type,eps_bc_low_type,eps_bc_high_type)
3272# eps
3273 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d=coeff(convw,convs,convl,vis3d,prand_eps,prand_eps,f1_sst,scheme_turb)
3274 su3d,sp3d,fmu3d= calceps(su3d,sp3d,k3d,eps3d,vis3d,gen)
3275
3276 ap3d,su3d=crank_nicol(eps3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv_keps)
3277
3278# this is needed only when blockage ios used
3279 if i_block_start !=0 or i_block_end !=0 or j_block_start !=0 or i_block_end !=0:
3280 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_block('eps')
3281
3282 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=fix_eps()
3283
3284 if solver_turb == 'pyamg':
3285 eps3d,residual_eps=solve_pyamg(eps3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_eps,acrank_conv_keps)
3286 elif solver_turb == 'pyamgx':
3287 eps3d,residual_eps=solve_pyamgx(eps3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_eps,acrank_conv_keps,'eps')
3288 else:
3289 eps3d,residual_eps=solve_3d(eps3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_eps,nsweep_keps,acrank_conv_keps,solver_turb,'eps')
3290
3291 eps3d=xp.maximum(eps3d,eps_min)
3292
3293 print(f"{'time eps: '}{time.time()-start_time:.2e}")
3294
3295 if pans or sst or kom_des or keps_des:
3296 fk3d=compute_fk(k3d,eps3d,fk3d)
3297
3298
3299# scale residuals
3300 residual_u=residual_u/resnorm_vel
3301 residual_v=residual_v/resnorm_vel
3302 residual_w=residual_w/resnorm_vel
3303 residual_p=residual_p/resnorm_p
3304 residual_k=residual_k/resnorm_vel**2
3305 residual_eps=residual_eps/resnorm_vel**3
3306 residual_om=residual_om/resnorm_vel
3307
3308# resmax=xp.max([residual_u ,residual_v,residual_p,residual_k,residual_eps,residual_om])
3309 resmax=xp.max(xp.array([residual_u ,residual_v,residual_w,residual_p]))
3310
3311 print(f"\n{'--time step:'}{itstep:4d}, {'iter:'}{iter:d}, {'max residual:'}{resmax:.2e}, {'u:'}{residual_u:.2e}\
3312, {'v:'}{residual_v:.2e}, {'w:'}{residual_w:.2e}, {'p:'}{residual_p:.2e}, {'k:'}{residual_k:.2e}\
3313, {'eps:'}{residual_eps:.2e}, {'om:'}{residual_om:.2e}\n")
3314
3315 print(f"\n{'monitor time step:'}{itstep:4d}, {'iter:'}{iter:1d}, {'u:'}{u3d[imon,jmon,kmon]: .2e}\
3316, {'v:'}{v3d[imon,jmon,kmon]: .2e}, {'w:'}{w3d[imon,jmon,kmon]: .2e}, {'p:'}{p3d[imon,jmon,kmon]: .2e}\
3317, {'k:'}{k3d[imon,jmon,kmon]: .2e}, {'eps:'}{eps3d[imon,jmon,kmon]: .2e}, {'om:'}{om3d[imon,jmon,kmon]: .2e}\n")
3318
3319
3320
3321 vismax=xp.max(vis3d.flatten())/viscos
3322 umax=xp.max(u3d.flatten())
3323 fkmax=xp.max(fk3d.flatten())
3324 epsmin=xp.min(eps3d.flatten())
3325 ommin=xp.min(om3d.flatten())
3326
3327
3328 kmin=xp.min(k3d.flatten())
3329
3330 if itstep%10 == 0:
3331 cfl_x=xp.abs(u3d)*dt[itstep]*areaw[1:,:,:]/vol
3332 cfl_y=xp.abs(v3d)*dt[itstep]*areas[:,1:,:]/vol
3333 cfl_x_max=xp.max(cfl_x)
3334 cfl_y_max=xp.max(cfl_y)
3335
3336# number of points larger than one
3337 cfl_y_no=cfl_y[xp.where( cfl_y > 1 )]
3338 cfl_x_no=cfl_x[xp.where( cfl_x > 1 )]
3339
3340 [i,j,k]= xp.where(cfl_y== xp.amax(cfl_y))
3341# check if i is an array
3342 if xp.size(i):
3343 ii=i[0]
3344 jj=j[0]
3345 else:
3346 ii=i
3347 jj=j
3348 print(f"\n{'-- cfl_max_x: '}{cfl_x_max:.2e}, {'cfl_max_y: '}{cfl_y_max:.2e}, {'at i= '}{ii:2d}, {'j= '}{jj:2d}\n")
3349
3350 print(f"\n{'No of points cfl_x > 1: '}{len(cfl_x_no):2d}, {'No of points cfl_y > 1: '}{len(cfl_y_no):2d}\n")
3351
3352
3353 print(f"\n{'--time step: '} {dt[itstep]:.2e}, {'iter: '}{iter:2d}, {'umax: '}{umax:.2e}, {'vismax: '}{vismax:.2e}, {'kmin: '}{kmin:.2e}, {'epsmin: '}{epsmin:.2e}, {'ommin: '}{ommin:.2e}, {'fkmax: '}{fkmax:.2e}\n")
3354
3355 print(f"{'time one iteration: '}{time.time()-start_time_iter:.2e}")
3356
3357# every 10th time step
3358 if itstep%10 == 0 and iter == 0:
3359 print(f"{'time every 10th time steps: '}{time.time()-time_10min:.2e}")
3360 time_10min = time.time()
3361
3362 if iter >= min_iter-1 and resmax < sormax:
3363
3364 break
3365
3366
3367 u3d_old,v3d_old,w3d_old,k3d_old,eps3d_old,om3d_old,dpdx_old,dpdy_old,dpdz_old=\
3368 update(u3d,v3d,w3d,k3d,eps3d,om3d,p3d_face_w,p3d_face_s,p3d_face_l)
3369# save data every itstep_save timsstep
3370 if itstep%itstep_save == 0 and itstep >= itstep_start:
3371 save_time_aver_data(u3d_mean,v3d_mean,w3d_mean,p3d_mean,eps3d_mean,om3d_mean,fk3d_mean,vis3d_mean,k3d_mean,uu3d_stress,\
3372 vv3d_stress,ww3d_stress,uv3d_stress,uw3d_stress,vw3d_stress,gen_mean)
3373 if save and itstep%itstep_save == 0 and itstep > 0:
3374 save_data(u3d,v3d,w3d,p3d,k3d,eps3d,om3d)
3375 if vtk and save_vtk_movie:
3376 save_vtk()
3377# save data if file save.data = 1
3378 if save and itstep%2 == 0:
3379 isave=0
3380 isave = xp.loadtxt('save.file')
3381 if isave == 1:
3382 save_data(u3d,v3d,w3d,p3d,k3d,eps3d,om3d)
3383 save_time_aver_data(u3d_mean,v3d_mean,w3d_mean,p3d_mean,eps3d_mean,om3d_mean,fk3d_mean,vis3d_mean,k3d_mean,uu3d_stress,\
3384 vv3d_stress,ww3d_stress,uv3d_stress,uw3d_stress,vw3d_stress,gen_mean)
3385
3386 if itstep >= itstep_start and itstep % itstep_stats == 0:
3387 u3d_mean,v3d_mean,w3d_mean,p3d_mean,k3d_mean,eps3d_mean,om3d_mean,uu3d_stress,vv3d_stress,ww3d_stress,uv3d_stress,uw3d_stress,vw3d_stress,\
3388 fk3d_mean,vis3d_mean,gen_mean= time_stats(u3d_mean,v3d_mean,w3d_mean,p3d_mean,k3d_mean,eps3d_mean,om3d_mean,uu3d_stress,\
3389 vv3d_stress,ww3d_stress,uv3d_stress,uw3d_stress,vw3d_stress,fk3d_mean,vis3d_mean,gen_mean)
3390
3391 vismean_mean=xp.max(vis3d_mean.flatten())/viscos/(itstep+1)
3392
3393 print('vismean_mean',vismean_mean)
3394
3395
3396# save data for restart
3397 if save:
3398 save_data(u3d,v3d,w3d,p3d,k3d,eps3d,om3d)
3399
3400 if vtk:
3401 itstep=ntstep
3402 save_vtk()
3403
3404# save time-averaged data
3405 save_time_aver_data(u3d_mean,v3d_mean,w3d_mean,p3d_mean,eps3d_mean,om3d_mean,fk3d_mean,vis3d_mean,k3d_mean,uu3d_stress,\
3406 vv3d_stress,ww3d_stress,uv3d_stress,uw3d_stress,vw3d_stress,gen_mean)
3407 print('program reached normal stop')
3408
3409main()
3410
solve_px(phi3d, aw3d, ae3d, as3d, an3d, al3d, ah3d, su3d, ap3d, tol_conv)
calck_kom(su3d, sp3d, k3d, om3d, vis3d, u3d_face_w, u3d_face_s, v3d_face_w, v3d_face_s, w3d_face_l)
compute_face_phi(phi3d, phi_bc_west, phi_bc_east, phi_bc_south, phi_bc_north, phi_bc_low, phi_bc_high, phi_bc_west_type, phi_bc_east_type, phi_bc_south_type, phi_bc_north_type, phi_bc_low_type, phi_bc_high_type, variable)
STG(U0, viscos, itstep, deltat, xp2d_synt, yp2d_synt, zp2d_synt, dw2d, dx2d, dy2d, dz2d, r11, r12, r13, r22, r23, r33, k_rans, om_rans)
read_restart_data(u3d, v3d, w3d, p3d, k3d, eps3d, om3d)
coeff_m(convw, convs, convl, vis3d, u3d, v3d, w3d)
vist_sst(vis3d, k3d, om3d, f2_sst, gen)
synt_fluct(nmodes, it, sli, yp, zp, uv_rans, visc, jmirror, dmin_synt)
solve_3d(phi3d, aw3d, ae3d, as3d, an3d, al3d, ah3d, su3d, ap3d, tol_conv, nmax, acrank_conv_local, solver_local, variable)
correct_conv(u3d, v3d, w3d, p3d, aw3d_p, as3d_p, al3d_p)
solve_pyamgx(phi3d, aw3d, ae3d, as3d, an3d, al3d, ah3d, su3d, ap3d, tol_conv, acrank_conv_local, variable)
vist_smag(u3d_face_w, u3d_face_s, u3d_face_l, v3d_face_w, v3d_face_s, v3d_face_l, w3d_face_w, w3d_face_s, w3d_face_l, vis3d)
calcu(su3d, sp3d, dpdx_old, p3d_face_w, p3d_face_s)
bc(su3d, sp3d, phi_bc_west, phi_bc_east, phi_bc_south, phi_bc_north, phi_bc_low, phi_bc_high, phi_bc_west_type, phi_bc_east_type, phi_bc_south_type, phi_bc_north_type, phi_bc_low_type, phi_bc_high_type)
calck(su3d, sp3d, k3d, eps3d, vis3d, u3d_face_w, u3d_face_s, v3d_face_w, v3d_face_s, w3d_face_l)
calcv(su3d, sp3d, dpdy_old, p3d_face_w, p3d_face_s)
calceps(su3d, sp3d, k3d, eps3d, vis3d, gen)
save_time_aver_data(u3d_mean, v3d_mean, w3d_mean, p3d_mean, eps3d_mean, om3d_mean, fk3d_mean, vis3d_mean, k3d_mean, uu3d_stress, vv3d_stress, ww3d_stress, uv3d_stress, uw3d_stress, vw3d_stress, gen_mean)
crank_nicol(phi3d_old, aw3d, ae3d, as3d, an3d, al3d, ah3d, apo3d, su3d, sp3d, acrank_conv_local)
calcp(convw, convs, convl)
save_data(u3d, v3d, w3d, p3d, k3d, eps3d, om3d)
dphidy(phi_face_w, phi_face_s)
update(u3d, v3d, w3d, k3d, eps3d, om3d, p3d_face_w, p3d_face_s, p3d_face_l)
dphidx(phi_face_w, phi_face_s)
calcom_sst(su3d, sp3d, om3d, gen, comm_term)
solve_p(phi3d, aw3d, ae3d, as3d, an3d, al3d, ah3d, su3d, ap3d, tol_conv)
solve_pyamg(phi3d, aw3d, ae3d, as3d, an3d, al3d, ah3d, su3d, ap3d, tol_conv, acrank_conv_local)
modify_conv(convw, convs, convl)
calcom(su3d, sp3d, om3d, gen, comm_term)
vist_wale(u3d_face_w, u3d_face_s, u3d_face_l, v3d_face_w, v3d_face_s, v3d_face_l, w3d_face_w, w3d_face_s, w3d_face_l, vis3d)
modify_init(u3d, v3d, w3d, k3d, om3d, eps3d, vis3d)
vist_pans(vis3d, k3d, eps3d, fmu3d)
modify_om(su3d, sp3d, comm_term)
coeff(convw, convs, convl, vis3d, prand_1, prand_2, f1_sst, scheme_local)
muscle_source(phi3d, cep, cem, cwp, cwm, cnp, cnm, csp, csm, chp, chm, clp, clm)
compute_fk(k3d, eps3d, fk3d)
time_stats(u3d_mean, v3d_mean, w3d_mean, p3d_mean, k3d_mean, eps3d_mean, om3d_mean, uu3d_stress, vv3d_stress, ww3d_stress, uv3d_stress, uw3d_stress, vw3d_stress, fk3d_mean, vis3d_mean, gen_mean)
modify_k(su3d, sp3d, gen)
vist_kom(vis3d, k3d, om3d)
conv(u3d, v3d, w3d, p3d_face_w, p3d_face_s, p3d_face_l)
calcw(su3d, sp3d, dpdz_old, p3d_face_l)