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
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
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
127 resnorm_p=uin*zmax*y2d[1,-1]
128 resnorm_vel=uin**2*zmax*y2d[1,-1]
138 jmirror_synt=int(nj/2)
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))
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))
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))
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))
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))
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))
227 eps_bc_south_type=
'd'
228 eps_bc_north_type=
'd'
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))
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
246 om_bc_south=xp.repeat(om_bc_south[:,
None], repeats=nk, axis=1)
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
254 om_bc_north=xp.repeat(om_bc_north[:,
None], repeats=nk, axis=1)
512def synt_fluct(nmodes,it,sli,yp,zp,uv_rans,visc,jmirror,dmin_synt):
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
530 uv_rans=xp.abs(uv_rans)
537 R=xp.genfromtxt(
"R.dat", dtype=
None,comments=
"%")
549 A=xp.genfromtxt(
"a.dat", dtype=
None,comments=
"%")
558 if xp.all(uv_rans==1):
561 up=(3.3*xp.max(uv_rans))**0.5
571 zp2d_synt=xp.repeat(zp[
None,:], repeats=nj, axis=0)
574 yp2d_synt=xp.repeat(yp[:,
None], repeats=nk, axis=1)
579 xp2d_synt=xp.ones((nj,nk))*xpp
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
588 dminy=xp.min(xp.diff(yp))
589 dminz=xp.min(xp.diff(zp))
590 dxmin=min(dminy,dminz)
593 dxmin=max(dxmin,dmin_synt)
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))
614 zp2d_wave=xp.zeros((nj,nk,nmodes+2))
621 wnrn=2.*math.pi/dxmin
624 wnre=9.*math.pi*amp/(55.*sli)
627 wnreta=(epsm/visc**3)**0.25
633 dkn=(wnrn-wnr1)/nmodes
636 wnr=xp.linspace(wnr1,wnrn,nmodes)
644 uv_rans_max=xp.max(uv_rans)
645 uv_rans_non=(uv_rans/(uv_rans_max+1e-10))**0.5
647 uv_rans_non=xp.repeat(uv_rans_non[:,
None], repeats=nk, axis=1)
650 print(f
"\n{'sli: '} {sli}, {'visc: '}{visc:.2e}, {'nmodes: '}{nmodes}, {'dxmin: '}{dxmin:.3e}, {'dkn: '}{dkn:.3e}, {'dmin_synt: '}{dmin_synt:.3e}")
652 print(f
"\n{'wnre: '} {wnre:.2e}, {'wnr1: '}{wnr1:.2e}, {'epsm: '}{epsm:.3e}, {'wnrn: '}{wnrn:.3e}")
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")
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)
672 print(
'time step no,',it)
676 kxio=xp.sin(teta)*xp.cos(fi)
677 kyio=xp.sin(teta)*xp.sin(fi)
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)
690 kxi=r11*kxio+r21*kyio+r31*kzio
691 kyi=r12*kxio+r22*kyio+r32*kzio
692 kzi=r13*kxio+r23*kyio+r33*kzio
694 sxi=r11*sxio+r21*syio+r31*szio
695 syi=r12*sxio+r22*syio+r32*szio
696 szi=r13*sxio+r23*syio+r33*szio
705 rk=xp.sqrt(kx**2+ky**2+kz**2)
709 xp2d_wave=xp.repeat(xp2d_synt[:,:,
None], repeats=nmodes, axis=2)
713 yp2d_wave=xp.repeat(yp2d_synt[:,:,
None], repeats=nmodes, axis=2)
716 zp2d_wave=xp.repeat(zp2d_synt[:,:,
None], repeats=nmodes, axis=2)
719 arg=arg1+arg2+arg3+psi
724 e=amp/wnre*(wnr/wnre)**4/((1.+(wnr/wnre)**2)**(17./6.))*xp.exp(-2*(wnr/wnreta)**2)
727 e=xp.where(rk < wnrn,e,0)
729 utn=xp.sqrt(e*up**2*dkn)
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)
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
742 uv=xp.mean(usynt_aniso*vsynt_aniso)
746 vsynt_aniso[jmirror:,:]=-vsynt_aniso[jmirror:,:]
749 uv_synt_mean=uv_synt_mean+uv
752 uvmean_time=xp.abs(uv_synt_mean)/(it+1)
753 print(
'uvmean_time',uvmean_time)
755 scale_2=(uv_rans_max/uvmean_time)**0.5*uv_rans_non
761 usynt_aniso=usynt_aniso*scale_2
762 vsynt_aniso=vsynt_aniso*scale_2
763 wsynt_aniso=wsynt_aniso*scale_2
767 uvmean_check=uvmean_check+xp.mean(usynt_aniso*vsynt_aniso,axis=1)
770 j=xp.where(uv_rans == xp.amax(uv_rans))
773 print(
'synt: uvmean',xp.abs(uvmean_check[j])/(it+1),
'uv_rans_max=',xp.max(xp.abs(uv_rans)),
'at j=',j)
775 return usynt_aniso,vsynt_aniso,wsynt_aniso
1339 if itstep == 0
and iter == 0:
1340 print(
'muscle scheme')
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)
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
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]
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,:,:]))
1365 visl[:,:,0]=0.5*(vis_turb[:,:,0]+vis_turb[:,:,-1])+viscos
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
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
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
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)
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)
1409 su3d_v=(1-blend)*su3d_v
1410 su3d_w=(1-blend)*su3d_w
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:,:,:]
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:,:]
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:]
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)))
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)))
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)))
1479 i1 = (xp.abs(x_embed-xp2d[:,0])).argmin()
1480 if itstep == 0
and iter == 0:
1481 print(
'embedded: i1',i1)
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:,:,:]
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:,:]
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:]
1511 if itstep == 0
and iter == 0:
1512 print(
'emnedded scheme; i_embed=',i1)
1515 apo3d=vol/dt[itstep]
1526 return aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,su3d_v,su3d_w
1584def coeff(convw,convs,convl,vis3d,prand_1,prand_2,f1_sst,scheme_local):
1586 if prand_1 == prand_2:
1587 prand = xp.ones((ni,nj,nk))*prand_1
1589 prand = f1_sst*prand_1+(1-f1_sst)*prand_2
1591 visw=xp.zeros((ni+1,nj,nk))
1592 viss=xp.zeros((ni,nj+1,nk))
1593 visl=xp.zeros((ni,nj,nk+1))
1595 vis_turb=(vis3d-viscos)/prand
1597 fk3d_local=xp.maximum(fk3d,fkmin_limit)
1598 vis_turb=(vis3d-viscos)/xp.abs(prand)/fk3d_local**2
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
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]
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,:,:]))
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]))
1621 if scheme_local ==
'h':
1622 if itstep == 0
and iter == 0:
1623 print(
'hybrid scheme, prand=',prand_1)
1625 aw3d=xp.maximum(convw[0:-1,:,:],diffw+(1-fx)*convw[0:-1,:,:])
1626 aw3d=xp.maximum(aw3d,0.)
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.)
1631 as3d=xp.maximum(convs[:,0:-1,:],diffs+(1-fy)*convs[:,0:-1,:])
1632 as3d=xp.maximum(as3d,0.)
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.)
1637 al3d=xp.maximum(convl[:,:,0:-1],diffl+(1-fz)*convl[:,:,0:-1])
1638 al3d=xp.maximum(al3d,0.)
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.)
1644 if scheme_local ==
'u':
1645 if itstep == 0
and iter == 0:
1646 print(
'upwind scheme, prand=',prand_1)
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)
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:,:,:]
1661 as3d=diffs+(1-fy)*convs[:,0:-1,:]
1662 an3d=xp.roll(diffs,-1,axis=1)-xp.roll(fy,-1,axis=1)*convs[:,1:,:]
1664 al3d=diffl+(1-fz)*convl[:,:,0:-1]
1665 ah3d=xp.roll(diffl,-1,axis=2)-xp.roll(fz,-1,axis=2)*convl[:,:,1:]
1668 apo3d=vol/dt[itstep]
1679 return aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d
2955 datax= xp.loadtxt(
"x2d.dat")
2958 datay= xp.loadtxt(
"y2d.dat")
2962 x2d=xp.zeros((ni+1,nj+1))
2963 y2d=xp.zeros((ni+1,nj+1))
2965 x2d=xp.reshape(x,(ni+1,nj+1))
2966 y2d=xp.reshape(y,(ni+1,nj+1))
2968 dataz=xp.loadtxt(
'z.dat')
2977 dz3d=xp.ones((ni+1,nj+1,nk))
2979 dz3d=xp.repeat(dz[
None,:], repeats=nj+1, axis=0)
2981 dz3d=xp.repeat(dz3d[
None,:,:], repeats=ni+1, axis=0)
2988 dz3d=xp.ones((ni+1,nj+1,nk))*dz
2989 z = xp.linspace(0, zmax, nk+1)
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:])
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))
3017 for itstep
in range(0,ntstep):
3020 for iter
in range(0,maxit):
3022 start_time_iter = time.time()
3024 start_time = time.time()
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()
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)
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
3035 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d=
coeff(convw,convs,convl,vis3d,1,1,f1_sst,scheme)
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)
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')
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')
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()
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)
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)
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')
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')
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}")
3074 start_time = time.time()
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)
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)
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')
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')
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}")
3096 start_time = time.time()
3098 convw,convs,convl=
conv(u3d,v3d,w3d,p3d_face_w,p3d_face_s,p3d_face_l)
3106 su3d=(convw[0:-1,:,:]-convw[1:,:,:]\
3107 +convs[:,0:-1,:]-convs[:,1:,:]\
3108 +convl[:,:,0:-1]-convl[:,:,1:])/acrank/dt[itstep]
3111 if iter == 0
and itstep == 0:
3112 aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d=
calcp(convw,convs,convl)
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')
3117 ap3d=aw3d+ae3d+as3d+an3d+al3d+ah3d
3123 print(
'xp.sum(su3d)',xp.sum(su3d))
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':
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')
3132 elif solver_p ==
'pyamgx_p':
3134 p3d=
solve_px(p3d,aw3d,ae3d,as3d,an3d,al3d,ah3d,su3d,ap3d,convergence_limit_p)
3136 print(
'solver_p = ',solver_p)
3137 print(
'no such solver')
3141 convw,convs,convl,p3d,u3d,v3d,w3d,su3d=
correct_conv(u3d,v3d,w3d,p3d,aw3d_p,as3d_p,al3d_p)
3143 res_1d=su3d.flatten()
3144 residual_p=xp.linalg.norm(res_1d,ord=1)
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')
3156 print(f
"{'time p: '}{time.time()-start_time:.2e}")
3158 start_time = time.time()
3160 if sst
or kom
or kom_des:
3162 vis3d=
vist_sst(vis3d,k3d,om3d,f2_sst,gen)
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)
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)
3172 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d=
coeff(convw,convs,convl,vis3d,prand_k,prand_k,f1_sst,scheme_turb)
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)
3177 ap3d,su3d=
crank_nicol(k3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv_kom)
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')
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')
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}")
3191 k3d=xp.maximum(k3d,k_min)
3193 start_time = time.time()
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)
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)
3201 aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d=
coeff(convw,convs,convl,vis3d,prand_omega,prand_omega,f1_sst,scheme_turb)
3203 su3d,sp3d,f1_sst,f2_sst=
calcom_sst(su3d,sp3d,om3d,gen,comm_term)
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)
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')
3212 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=
fix_omega()
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')
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)
3222 print(f
"{'time omega: '}{time.time()-start_time:.2e}")
3224 start_time = time.time()
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)
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)
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:
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)
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)
3247 ap3d,su3d=
crank_nicol(k3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv_keps)
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')
3253 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=
fix_k()
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')
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')
3263 k3d=xp.maximum(k3d,k_min)
3264 print(f
"{'time k: '}{time.time()-start_time:.2e}")
3265 start_time = time.time()
3268 start_time = time.time()
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)
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)
3276 ap3d,su3d=
crank_nicol(eps3d_old,aw3d,ae3d,as3d,an3d,al3d,ah3d,apo3d,su3d,sp3d,acrank_conv_keps)
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')
3282 aw3d,ae3d,as3d,an3d,al3d,ah3d,ap3d,su3d,sp3d=
fix_eps()
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')
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')
3291 eps3d=xp.maximum(eps3d,eps_min)
3293 print(f
"{'time eps: '}{time.time()-start_time:.2e}")
3295 if pans
or sst
or kom_des
or keps_des:
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
3309 resmax=xp.max(xp.array([residual_u ,residual_v,residual_w,residual_p]))
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")
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")
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())
3328 kmin=xp.min(k3d.flatten())
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)
3337 cfl_y_no=cfl_y[xp.where( cfl_y > 1 )]
3338 cfl_x_no=cfl_x[xp.where( cfl_x > 1 )]
3340 [i,j,k]= xp.where(cfl_y== xp.amax(cfl_y))
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")
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")
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")
3355 print(f
"{'time one iteration: '}{time.time()-start_time_iter:.2e}")
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()
3362 if iter >= min_iter-1
and resmax < sormax:
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)
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:
3378 if save
and itstep%2 == 0:
3380 isave = xp.loadtxt(
'save.file')
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)
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)
3391 vismean_mean=xp.max(vis3d_mean.flatten())/viscos/(itstep+1)
3393 print(
'vismean_mean',vismean_mean)
3398 save_data(u3d,v3d,w3d,p3d,k3d,eps3d,om3d)
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')