8 global c_omega_1, c_omega_2, cmu, convergence_limit_eps, convergence_limit_k, convergence_limit_om, convergence_limit_pp, \
9 convergence_limit_u, convergence_limit_v, convergence_limit_w, cyclic_x, dist, earsm,fx, fy,imon,jmon,kappa,k_bc_east,k_bc_east_type, \
10 k_bc_north,k_bc_north_type,k_bc_south, k_bc_south_type,k_bc_west,k_bc_west_type,kom,maxit, \
11 ni,nj,nsweep_kom, nsweep_pp, nsweep_vel, om_bc_east, om_bc_east_type, om_bc_north, om_bc_north_type, \
12 om_bc_south, om_bc_south_type, om_bc_west, om_bc_west_type, p_bc_east, p_bc_east_type, \
13 p_bc_north, p_bc_north_type, p_bc_south, p_bc_south_type, p_bc_west, p_bc_west_type, pinn, \
14 prand_k,prand_omega,resnorm_p,resnorm_vel,restart,save,save_vtk_movie,scheme,scheme_turb,solver_pp,solver_vel, \
15 solver_turb,sormax,u_bc_east,u_bc_east_type,u_bc_north,u_bc_north_type,u_bc_south,u_bc_south_type,u_bc_west, \
16 u_bc_west_type,urfvis,urf_vel,urf_k,urf_p,urf_omega,v_bc_east,v_bc_east_type,v_bc_north,v_bc_north_type, \
17 v_bc_south,v_bc_south_type,v_bc_west,v_bc_west_type,viscos,vol,tk,tk_save,tk_file_name,\
18 uu_bc_west,uu_bc_east,uu_bc_south,uu_bc_north,uu_bc_west_type,uu_bc_east_type,uu_bc_south_type,uu_bc_north_type,\
19 uv_bc_west,uv_bc_east,uv_bc_south,uv_bc_north,uv_bc_west_type,uv_bc_east_type,uv_bc_south_type,uv_bc_north_type,\
20 vv_bc_west,vv_bc_east,vv_bc_south,vv_bc_north,vv_bc_west_type,vv_bc_east_type,vv_bc_south_type,vv_bc_north_type, \
70 convergence_limit_u=-1e-6
71 convergence_limit_v=1e-6
72 convergence_limit_k=-1e-10
73 convergence_limit_om=-1e-10
74 convergence_limit_pp=5e-4
88 resnorm_p=uin*y2d[1,-1]
89 resnorm_vel=uin**2*y2d[1,-1]
96 u_bc_east=np.zeros(nj)
97 u_bc_south=np.zeros(ni)
98 u_bc_north=np.zeros(ni)
106 v_bc_west=np.zeros(nj)
107 v_bc_east=np.zeros(nj)
108 v_bc_south=np.zeros(ni)
109 v_bc_north=np.zeros(ni)
117 p_bc_west=np.zeros(nj)
118 p_bc_east=np.zeros(nj)
119 p_bc_south=np.zeros(ni)
120 p_bc_north=np.zeros(ni)
128 k_bc_west=np.zeros(nj)
129 k_bc_east=np.zeros(nj)
130 k_bc_south=np.zeros(ni)
131 k_bc_north=np.zeros(ni)
139 om_bc_west=np.zeros(nj)
140 om_bc_east=np.zeros(nj)
142 xwall_s=0.5*(x2d[0:-1,0]+x2d[1:,0])
143 ywall_s=0.5*(y2d[0:-1,0]+y2d[1:,0])
144 dist2_s=(yp2d[:,0]-ywall_s)**2+(xp2d[:,0]-xwall_s)**2
145 om_bc_south=6*viscos/0.075/dist2_s
147 xwall_n=0.5*(x2d[0:-1,-1]+x2d[1:,-1])
148 ywall_n=0.5*(y2d[0:-1,-1]+y2d[1:,-1])
149 dist2_n=(yp2d[:,-1]-ywall_n)**2+(xp2d[:,-1]-xwall_n)**2
150 om_bc_north=6*viscos/0.075/dist2_n
726def solve_2d(phi2d,aw2d,ae2d,as2d,an2d,su2d,ap2d,tol_conv,nmax,solver_local):
728 print(
'solve_2d called')
731 aw=np.matrix.flatten(aw2d)
732 ae=np.matrix.flatten(ae2d)
733 as1=np.matrix.flatten(as2d)
734 an=np.matrix.flatten(an2d)
735 ap=np.matrix.flatten(ap2d)
742 A = sparse.diags([ap, -an[:-1], -as1[1:], -ae, -aw[nj:],-aw,-ae[nj*(ni-1):]],\
743 [0, 1, -1, nj, -nj,nj*(ni-1), -nj*(ni-1)], format=
'csr')
745 A = sparse.diags([ap, -an[0:-1], -as1[1:], -ae, -aw[nj:]], [0, 1, -1, nj, -nj], format=
'csr')
747 su=np.matrix.flatten(su2d)
748 phi=np.matrix.flatten(phi2d)
750 res_su=np.linalg.norm(su)
751 resid_init=np.linalg.norm(A*phi - su)
765 abs_tol =abs(tol_conv)*resid_init
768 if solver_local ==
'direct':
770 print(
'solver in solve_2d: direct solver')
772 resid=np.linalg.norm(A*phi - su)
773 phi = linalg.spsolve(A,su)
774 if solver_local ==
'pyamg':
776 print(
'solver in solve_2d: pyamg solver')
777 App = pyamg.ruge_stuben_solver(A)
779 phi = App.solve(su, tol=tol_conv, x0=phi, residuals=res_amg)
781 print(
'Residual history in pyAMG', [
"%0.4e" % i
for i
in res_amg])
782 if solver_local ==
'cgs':
784 print(
'solver in solve_2d: cgs')
785 phi,info=linalg.cgs(A,su,x0=phi, rtol=tol_conv, atol=abs_tol, maxiter=nmax)
786 if solver_local ==
'cg':
788 print(
'solver in solve_2d: cg')
789 phi,info=linalg.cg(A,su,x0=phi, rtol=tol_conv, atol=abs_tol, maxiter=nmax)
790 if solver_local ==
'gmres':
792 print(
'solver in solve_2d: gmres')
793 phi,info=linalg.gmres(A,su,x0=phi, rtol=tol_conv, atol=abs_tol, maxiter=nmax)
794 if solver_local ==
'qmr':
796 print(
'solver in solve_2d: qmr')
797 phi,info=linalg.qmr(A,su,x0=phi, rtol=tol_conv, atol=abs_tol, maxiter=nmax)
798 if solver_local ==
'lgmres':
800 print(
'solver in solve_2d: lgmres')
801 phi,info=linalg.lgmres(A,su,x0=phi, rtol=tol_conv, atol=abs_tol, maxiter=nmax)
803 print(
'warning in module solve_2d: convergence in sparse matrix solver not reached')
806 if solver_local !=
'direct':
807 resid=np.linalg.norm(A*phi - su)
809 delta_phi=np.max(np.abs(phi-phi_org))
811 phi2d=np.reshape(phi,(ni,nj))
812 phi2d_org=np.reshape(phi_org,(ni,nj))
814 if solver_local !=
'pyamg':
815 print(f
"{'residual history in solve_2d: initial residual: '} {resid_init:.2e}{'final residual: ':>30}{resid:.2e}\
816 {'delta_phi: ':>25}{delta_phi:.2e}")
819 return phi2d,resid_init
1148 datax= np.loadtxt(
"x2d.dat")
1151 datay= np.loadtxt(
"y2d.dat")
1155 x2d=np.zeros((ni+1,nj+1))
1156 y2d=np.zeros((ni+1,nj+1))
1158 x2d=np.reshape(x,(ni+1,nj+1))
1159 y2d=np.reshape(y,(ni+1,nj+1))
1162 xp2d=0.25*(x2d[0:-1,0:-1]+x2d[0:-1,1:]+x2d[1:,0:-1]+x2d[1:,1:])
1163 yp2d=0.25*(y2d[0:-1,0:-1]+y2d[0:-1,1:]+y2d[1:,0:-1]+y2d[1:,1:])
1167 vol=np.zeros((ni,nj))
1168 areas=np.zeros((ni,nj+1))
1169 areasx=np.zeros((ni,nj+1))
1170 areasy=np.zeros((ni,nj+1))
1171 areaw=np.zeros((ni+1,nj))
1172 areawx=np.zeros((ni+1,nj))
1173 areawy=np.zeros((ni+1,nj))
1174 areaz=np.zeros((ni,nj))
1175 as_bound=np.zeros((ni))
1176 an_bound=np.zeros((ni))
1177 aw_bound=np.zeros((nj))
1178 ae_bound=np.zeros((nj))
1179 az_bound=np.zeros((ni,nj))
1180 fx=np.zeros((ni,nj))
1181 fy=np.zeros((ni,nj))
1185 u_bc_west=np.ones(nj)
1186 u_bc_east=np.zeros(nj)
1187 u_bc_south=np.zeros(ni)
1188 u_bc_north=np.zeros(ni)
1196 v_bc_west=np.zeros(nj)
1197 v_bc_east=np.zeros(nj)
1198 v_bc_south=np.zeros(ni)
1199 v_bc_north=np.zeros(ni)
1207 p_bc_west=np.zeros(nj)
1208 p_bc_east=np.zeros(nj)
1209 p_bc_south=np.zeros(ni)
1210 p_bc_north=np.zeros(ni)
1218 k_bc_west=np.zeros(nj)
1219 k_bc_east=np.zeros(nj)
1220 k_bc_south=np.zeros(ni)
1221 k_bc_north=np.zeros(ni)
1229 om_bc_west=np.zeros(nj)
1230 om_bc_east=np.zeros(nj)
1231 om_bc_south=np.zeros(ni)
1232 om_bc_north=np.zeros(ni)
1236 om_bc_south_type=
'd'
1237 om_bc_north_type=
'd'
1240 uu_bc_west=np.zeros(nj)
1241 uu_bc_east=np.zeros(nj)
1242 uu_bc_south=np.zeros(ni)
1243 uu_bc_north=np.zeros(ni)
1247 uu_bc_south_type=
'd'
1248 uu_bc_north_type=
'd'
1251 uv_bc_west=np.zeros(nj)
1252 uv_bc_east=np.zeros(nj)
1253 uv_bc_south=np.zeros(ni)
1254 uv_bc_north=np.zeros(ni)
1258 uv_bc_south_type=
'd'
1259 uv_bc_north_type=
'd'
1262 vv_bc_west=np.zeros(nj)
1263 vv_bc_east=np.zeros(nj)
1264 vv_bc_south=np.zeros(ni)
1265 vv_bc_north=np.zeros(ni)
1269 vv_bc_south_type=
'd'
1270 vv_bc_north_type=
'd'
1285 areaw,areawx,areawy,areas,areasx,areasy,vol,fx,fy,aw_bound,ae_bound,as_bound,an_bound,dist=
init()
1289 u2d=np.ones((ni,nj))*1e-20
1290 v2d=np.ones((ni,nj))*1e-20
1291 p2d=np.ones((ni,nj))*1e-20
1292 pp2d=np.ones((ni,nj))*1e-20
1293 k2d=np.ones((ni,nj))*1
1294 om2d=np.ones((ni,nj))*1
1295 vis2d=np.ones((ni,nj))*viscos
1296 vis2d_earsm=np.ones((ni,nj))*viscos
1297 uu2d=np.ones((ni,nj))*1e-20
1298 uv2d=np.ones((ni,nj))*1e-20
1299 vv2d=np.ones((ni,nj))*1e-20
1300 ww2d=np.ones((ni,nj))*1e-20
1302 fmu2d=np.ones((ni,nj))
1303 gen=np.ones((ni,nj))
1305 convw=np.ones((ni+1,nj))*1e-20
1306 convs=np.ones((ni,nj+1))*1e-20
1308 aw2d=np.ones((ni,nj))*1e-20
1309 ae2d=np.ones((ni,nj))*1e-20
1310 as2d=np.ones((ni,nj))*1e-20
1311 an2d=np.ones((ni,nj))*1e-20
1312 al2d=np.ones((ni,nj))*1e-20
1313 ah2d=np.ones((ni,nj))*1e-20
1314 ap2d=np.ones((ni,nj))*1e-20
1315 ap2d_vel=np.ones((ni,nj))*1e-20
1316 su2d=np.ones((ni,nj))*1e-20
1317 sp2d=np.ones((ni,nj))*1e-20
1318 ap2d=np.ones((ni,nj))*1e-20
1319 dudx=np.ones((ni,nj))*1e-20
1320 dudy=np.ones((ni,nj))*1e-20
1321 usynt_inlet=np.ones((nj))*1e-20
1322 vsynt_inlet=np.ones((nj))*1e-20
1323 wsynt_inlet=np.ones((nj))*1e-20
1325 uu2d_face_w = np.ones((ni,nj))*1e-20
1326 uu2d_face_s = np.ones((ni,nj))*1e-20
1327 vv2d_face_w = np.ones((ni,nj))*1e-20
1328 vv2d_face_s = np.ones((ni,nj))*1e-20
1329 uv2d_face_w = np.ones((ni,nj))*1e-20
1330 uv2d_face_s = np.ones((ni,nj))*1e-20
1333 delta_max=np.maximum(vol/areas[:,1:],vol/areaw[1:,:])
1338 c_k =np.ones((ni,nj))
1339 prand_k =np.ones((ni,nj))*2
1340 c_omega_2 =np.ones((ni,nj))*c_omega_2
1349 u2d,v2d,k2d,om2d,vis2d,dist=
modify_init(u2d,v2d,k2d,om2d,vis2d)
1355 k2d=np.maximum(k2d,1e-6)
1357 u2d_face_w,u2d_face_s=
compute_face_phi(u2d,u_bc_west,u_bc_east,u_bc_south,u_bc_north,\
1358 u_bc_west_type,u_bc_east_type,u_bc_south_type,u_bc_north_type)
1359 v2d_face_w,v2d_face_s=
compute_face_phi(v2d,v_bc_west,v_bc_east,v_bc_south,v_bc_north,\
1360 v_bc_west_type,v_bc_east_type,v_bc_south_type,v_bc_north_type)
1361 p2d_face_w,p2d_face_s=
compute_face_phi(p2d,p_bc_west,p_bc_east,p_bc_south,p_bc_north,\
1362 p_bc_west_type,p_bc_east_type,p_bc_south_type,p_bc_north_type)
1366 u_bc_west,v_bc_west,k_bc_west,om_bc_west,u2d_face_w,convw =
modify_inlet()
1368 convw,convs=
conv(u2d,v2d,p2d_face_w,p2d_face_s)
1389 for iter
in range(0,abs(maxit)):
1391 start_time_iter = time.time()
1393 start_time = time.time()
1397 u_bc_west,v_bc_west,k_bc_west,om_bc_west,u2d_face_w,convw =
modify_inlet()
1399 uu2d,vv2d,ww2d,uv2d,vis2d_earsm=calc_earsm(k2d,om2d,u2d_face_w,u2d_face_s,v2d_face_w,v2d_face_s,\
1400 uu2d,uv2d,vv2d,ww2d,vis2d_earsm)
1401 aw2d,ae2d,as2d,an2d,su2d,sp2d=
coeff(convw,convs,vis2d_earsm,1,scheme)
1403 aw2d,ae2d,as2d,an2d,su2d,sp2d=
coeff(convw,convs,vis2d,1,scheme)
1406 uu2d_face_w,uu2d_face_s=
compute_face_phi(uu2d,uu_bc_west,uu_bc_east,uu_bc_south,uu_bc_north,\
1407 uu_bc_west_type,uu_bc_east_type,uu_bc_south_type,uu_bc_north_type)
1408 uv2d_face_w,uv2d_face_s=
compute_face_phi(uv2d,uv_bc_west,uv_bc_east,uv_bc_south,uv_bc_north,\
1409 uv_bc_west_type,uv_bc_east_type,uv_bc_south_type,uv_bc_north_type)
1410 vv2d_face_w,vv2d_face_s=
compute_face_phi(vv2d,vv_bc_west,vv_bc_east,vv_bc_south,vv_bc_north,\
1411 vv_bc_west_type,vv_bc_east_type,vv_bc_south_type,vv_bc_north_type)
1415 su2d,sp2d=
bc(su2d,sp2d,u_bc_west,u_bc_east,u_bc_south,u_bc_north, \
1416 u_bc_west_type,u_bc_east_type,u_bc_south_type,u_bc_north_type)
1417 su2d,sp2d,ap2d=
calcu(su2d,sp2d,p2d_face_w,p2d_face_s,uu2d_face_w,uu2d_face_s,uv2d_face_w,uv2d_face_s)
1420 u2d,residual_u=
solve_2d(u2d,aw2d,ae2d,as2d,an2d,su2d,ap2d,convergence_limit_u,nsweep_vel,solver_vel)
1421 print(f
"{'time u: '}{time.time()-start_time:.2e}")
1423 start_time = time.time()
1426 su2d,sp2d=
bc(su2d,sp2d,v_bc_west,v_bc_east,v_bc_south,v_bc_north, \
1427 v_bc_west_type,v_bc_east_type,v_bc_south_type,v_bc_north_type)
1428 su2d,sp2d,ap2d,ap2d_vel=
calcv(su2d,sp2d,p2d_face_w,p2d_face_s,uv2d_face_w,uv2d_face_s,vv2d_face_w,vv2d_face_s)
1430 v2d,residual_v=
solve_2d(v2d,aw2d,ae2d,as2d,an2d,su2d,ap2d,convergence_limit_v,nsweep_vel,solver_vel)
1431 print(f
"{'time v: '}{time.time()-start_time:.2e}")
1433 start_time = time.time()
1435 convw,convs=
conv(u2d,v2d,p2d_face_w,p2d_face_s)
1438 aw2d,ae2d,as2d,an2d,su2d,ap2d=
calcp(pp2d,ap2d_vel)
1439 pp2d=np.zeros((ni,nj))
1441 pp2d,residual_p=
solve_2d(pp2d,aw2d,ae2d,as2d,an2d,su2d,ap2d,convergence_limit_pp,nsweep_pp,solver_pp)
1448 su2d=convw[0:-1,:]-np.roll(convw[0:-1,:],-1,axis=0)+convs[:,0:-1]-np.roll(convs[:,0:-1],-1,axis=1)
1449 residual_pp=abs(np.sum(su2d))
1451 print(f
"{'time pp: '}{time.time()-start_time:.2e}")
1453 u2d_face_w,u2d_face_s=
compute_face_phi(u2d,u_bc_west,u_bc_east,u_bc_south,u_bc_north,\
1454 u_bc_west_type,u_bc_east_type,u_bc_south_type,u_bc_north_type)
1455 v2d_face_w,v2d_face_s=
compute_face_phi(v2d,v_bc_west,v_bc_east,v_bc_south,v_bc_north,\
1456 v_bc_west_type,v_bc_east_type,v_bc_south_type,v_bc_north_type)
1457 p2d_face_w,p2d_face_s=
compute_face_phi(p2d,p_bc_west,p_bc_east,p_bc_south,p_bc_north,\
1458 p_bc_west_type,p_bc_east_type,p_bc_south_type,p_bc_north_type)
1460 start_time = time.time()
1467 start_time = time.time()
1469 aw2d,ae2d,as2d,an2d,su2d,sp2d=
coeff(convw,convs,vis2d_earsm,prand_k,scheme_turb)
1471 aw2d,ae2d,as2d,an2d,su2d,sp2d=
coeff(convw,convs,vis2d,prand_k,scheme_turb)
1474 su2d,sp2d=
bc(su2d,sp2d,k_bc_west,k_bc_east,k_bc_south,k_bc_north, \
1475 k_bc_west_type,k_bc_east_type,k_bc_south_type,k_bc_north_type)
1476 su2d,sp2d,gen,ap2d=
calck(su2d,sp2d,k2d,om2d,vis2d,u2d_face_w,u2d_face_s,v2d_face_w,v2d_face_s)
1478 aw3d,ae3d,as3d,an3d,ap3d,su3d,sp3d=
fix_k()
1481 k2d,residual_k=
solve_2d(k2d,aw2d,ae2d,as2d,an2d,su2d,ap2d,convergence_limit_k,nsweep_kom,solver_turb)
1482 k2d=np.maximum(k2d,1e-10)
1483 print(f
"{'time k: '}{time.time()-start_time:.2e}")
1485 start_time = time.time()
1489 aw2d,ae2d,as2d,an2d,su2d,sp2d=
coeff(convw,convs,vis2d_earsm,prand_omega,scheme_turb)
1491 aw2d,ae2d,as2d,an2d,su2d,sp2d=
coeff(convw,convs,vis2d,prand_omega,scheme_turb)
1492 su2d,sp2d=
bc(su2d,sp2d,om_bc_west,om_bc_east,om_bc_south,om_bc_north,\
1493 om_bc_west_type,om_bc_east_type,om_bc_south_type,om_bc_north_type)
1494 su2d,sp2d,ap2d=
calcom(su2d,sp2d,om2d,gen)
1496 aw2d,ae2d,as2d,an2d,ap2d,su2d,sp2d=
fix_omega()
1499 om2d,residual_om=
solve_2d(om2d,aw2d,ae2d,as2d,an2d,su2d,ap2d,convergence_limit_om,nsweep_kom,solver_turb)
1500 om2d=np.maximum(om2d,1e-10)
1502 print(f
"{'time omega: '}{time.time()-start_time:.2e}")
1505 residual_u=residual_u/resnorm_vel
1506 residual_v=residual_v/resnorm_vel
1507 residual_pp=residual_p/resnorm_p
1508 residual_k=residual_k/resnorm_vel**2
1509 residual_om=residual_om/resnorm_vel
1511 resmax=np.max([residual_u ,residual_v,residual_pp])
1513 print(f
"\n{'--iter:'}{iter:d}, {'max residual:'}{resmax:.2e}, {'u:'}{residual_u:.2e}\
1514, {'v:'}{residual_v:.2e}, {'cont:'}{residual_pp:.2e}, {'k:'}{residual_k:.2e}\
1515, {'om:'}{residual_om:.2e}\n")
1517 print(f
"\n{'monitor iteration:'}{iter:4d}, {'u:'}{u2d[imon,jmon]: .2e}\
1518, {'v:'}{v2d[imon,jmon]: .2e}, {'p:'}{p2d[imon,jmon]: .2e}\
1519, {'k:'}{k2d[imon,jmon]: .2e}, {'om:'}{om2d[imon,jmon]: .2e}\n")
1523 vismax=np.max(vis2d.flatten())/viscos
1524 umax=np.max(u2d.flatten())
1525 ommin=np.min(om2d.flatten())
1526 kmin=np.min(k2d.flatten())
1528 print(f
"\n{'---iter: '}{iter:2d}, {'umax: '}{umax:.2e},{'vismax: '}{vismax:.2e}, {'kmin: '}{kmin:.2e}, {'ommin: '}{ommin:.2e}\n")
1530 print(f
"{'time one iteration: '}{time.time()-start_time_iter:.2e}")
1532 if resmax < sormax
and iter > 0:
1539 save_data(u2d,v2d,p2d,k2d,om2d,vis2d,ap2d_vel)
1545 print(
'program reached normal stop')