pyCALC-RANS
Loading...
Searching...
No Matches
exec-pyCALC-RANS.py
Go to the documentation of this file.
1from scipy import sparse
2import numpy as np
3import sys
4import time
5from scipy.sparse import spdiags,linalg,eye
6
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, \
21 x2d,xp2d,y2d,yp2d
22
23
24
25
26
27 import numpy as np
28 import sys
29
30
31# section 1 choice of differencing scheme #
32 scheme='h' #hybrid
33 scheme_turb='h' #hybrid upwind-central
34
35# section 2 turbulence models #
36 cmu=0.09
37 kom = True
38 c_omega_1= 5./9.
39 c_omega_2=3./40.
40 prand_omega=2.0
41 prand_k=2.0
42 pinn = True
43
44# section 3 restart/save #
45 restart = True
46 save = True
47
48# section 4 fluid properties #
49 viscos=1/2000
50
51# section 5 relaxation factors #
52 urfvis=0.5
53 urf_vel=0.5
54 urf_k=0.5
55 urf_p=1.0
56 urf_omega=0.5
57
58# section 6 number of iteration and convergence criterira #
59 maxit=20000
60 min_iter=1
61 sormax=1e-8
62
63 solver_vel='lgmres'
64 solver_pp='direct'
65 solver_turb='direct'
66
67 nsweep_vel=50
68 nsweep_pp=50
69 nsweep_kom=50
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
75
76
77# section 7 all variables are printed during the iteration at node #
78 imon=0
79 jmon=30
80
81# section 8 save data for post-processing #
82 vtk=False
83 save_all_files=False
84 vtk_file_name='bound'
85
86# section 9 residual scaling parameters #
87 uin=20
88 resnorm_p=uin*y2d[1,-1]
89 resnorm_vel=uin**2*y2d[1,-1]
90
91
92# Section 10 boundary conditions #
93
94# boundary conditions for u
95 u_bc_west=np.ones(nj)
96 u_bc_east=np.zeros(nj)
97 u_bc_south=np.zeros(ni)
98 u_bc_north=np.zeros(ni)
99
100 u_bc_west_type='n'
101 u_bc_east_type='n'
102 u_bc_south_type='d'
103 u_bc_north_type='n'
104
105# boundary conditions for v
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)
110
111 v_bc_west_type='n'
112 v_bc_east_type='n'
113 v_bc_south_type='d'
114 v_bc_north_type='d'
115
116# boundary conditions for p
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)
121
122 p_bc_west_type='n'
123 p_bc_east_type='n'
124 p_bc_south_type='n'
125 p_bc_north_type='n'
126
127# boundary conditions for k
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)
132
133 k_bc_west_type='n'
134 k_bc_east_type='n'
135 k_bc_south_type='d'
136 k_bc_north_type='n'
137
138# boundary conditions for omega
139 om_bc_west=np.zeros(nj)
140 om_bc_east=np.zeros(nj)
141
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
146
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
151
152 om_bc_west_type='n'
153 om_bc_east_type='n'
154 om_bc_south_type='d'
155 om_bc_north_type='n'
156
157 return
158
159
160def modify_init(u2d,v2d,k2d,om2d,vis2d):
161
162# set inlet field in entre domain
163 u2d=np.repeat(u_bc_west[None,:], repeats=ni, axis=0)
164 k2d=np.ones((ni,nj))
165 om2d=np.ones((ni,nj))
166
167 vis2d=k2d/om2d+viscos
168
169 return u2d,v2d,k2d,om2d,vis2d,dist
170
172
173# load vist/viscos predicted by PINN
174 vist_ML_5200 = np.loadtxt('../PINN/vist_pred-PINN-from-vist-diffusion-pinn-5200-plus-units-load-5-cells-fixed.txt')
175
176 viscos_5200 = 1/5200
177
178 vist_ML_5200 = vist_ML_5200*viscos_5200
179
180 data_5200= np.loadtxt('y_u_k_om_uv_5200-RANS-half-channel.txt')
181 y_5200 = data_5200[:,0]
182 y1d= yp2d[0,:]
183
184 prand_k_ML_5200 = np.loadtxt('../PINN/prand_k_5200-plus-units-from-balance.txt')
185
186 print('prand_k_ML_5200',prand_k_ML_5200)
187
188 vist_ML = np.interp(y1d, y_5200, vist_ML_5200)
189 prand_k = np.interp(y1d, y_5200, prand_k_ML_5200)
190
191# make it 2D
192 prand_k = np.repeat(prand_k[None,:], repeats=ni, axis=0)
193
194 print('vist_ML',vist_ML)
195
196 vist = k2d/om2d
197
198 print('prand_k',prand_k)
199
200# load c_k taken from balance of k eq.
201 c_k_ML_5200= np.loadtxt('../PINN/c_k_pred_5200-plus-units-from-balance.txt')
202 c_k_ML = np.interp(y1d, y_5200, c_k_ML_5200)
203# make it 2D
204 c_k = np.repeat(c_k_ML[None,:], repeats=ni, axis=0)
205
206# load c_omega_2 taken from balance of omega eq.
207 c_omega_2_ML_5200 = np.loadtxt('../PINN/c_omega_2_pred_5200-plus-units-from-balance.txt')
208 c_omega_2_ML = np.interp(y1d, y_5200, c_omega_2_ML_5200)
209# make it 2D
210 c_omega_2 = np.repeat(c_omega_2_ML[None,:], repeats=ni, axis=0)
211
212 return prand_k, c_k, c_omega_2
213
215
216 global y_rans,y_rans,u_rans,v_rans,k_rans,om_rans,uv_rans,k_bc_west,eps_bc_west,om_bc_west
217
218 return u_bc_west,v_bc_west,k_bc_west,om_bc_west,u2d_face_w,convw
219
220def modify_conv(convw,convs):
221
222# since we are solving for fully-developed channel flow, we know that the convection terms are zero
223 convs=np.zeros((ni,nj))
224 convw=np.zeros((ni,nj))
225
226 return convw,convs
227
228def modify_u(su2d,sp2d):
229
230 global file1
231
232# add a driving pressure gradient term
233 su2d= su2d+vol
234
235# we know that the convection and diffusion term in the x direction are zero
236 aw2d=np.zeros((ni,nj))
237 ae2d=np.zeros((ni,nj))
238
239# we know that for this flow the wall shear stress mustt be equal to one (since the driving pressure
240# gradient is equal to one). We print it every iteration to see if it is one. When it reaches one it is
241# a good indicator that the flow has converged
242
243 tauw_south=viscos*np.sum(as_bound*u2d[:,0])/x2d[-1,0]
244 tauw_north=viscos*np.sum(an_bound*u2d[:,-1])/x2d[-1,0]
245
246 print(f"{'tau wall, south: '} {tauw_south:.3f},{' tau wall, north: '} {tauw_north:.3f}")
247
248
249 return su2d,sp2d
250
251def modify_v(su2d,sp2d):
252
253 return su2d,sp2d
254
255def modify_p(su2d,sp2d):
256
257 return su2d,sp2d
258
259def modify_k(su2d,sp2d):
260
261 if iter == 0:
262 np.savetxt('k-iteration.dat', np.c_[iter,k2d[-1,5],k2d[-1,10],k2d[-1,20],k2d[-1,30],k2d[-1,40],\
263 k2d[-1,50],k2d[-1,58]])
264 else:
265 with open('k-iteration.dat','ab') as f:
266 np.savetxt(f,np.c_[iter,k2d[-1,5],k2d[-1,10],k2d[-1,20],k2d[-1,30],k2d[-1,40],\
267 k2d[-1,50],k2d[-1,58]])
268
269 return su2d,sp2d
270
271def modify_om(su2d,sp2d):
272
273 return su2d,sp2d
274
275def modify_outlet(convw):
276
277# since we are solving for fully-developed channel flow, we know that the convection terms are zero
278 convw=np.zeros((ni+1,nj))
279
280 return convw
281
283
284 aw2d[:,0]=0
285 ae2d[:,0]=0
286 as2d[:,0]=0
287 an2d[:,0]=0
288 ap2d[:,0]=1
289 su2d[:,0]=om_bc_south
290
291 return aw2d,ae2d,as2d,an2d,ap2d,su2d,sp2d
292
293def modify_vis(vis2d):
294
295 return vis2d
296
297
298def fix_k():
299
300 return aw2d,ae2d,as2d,an2d,ap2d,su2d,sp2d
301from scipy import sparse
302import numpy as np
303import sys
304import time
305import pyamg
306from scipy.sparse import spdiags,linalg,eye
307import socket
308
309def init():
310 print('hostname: ',socket.gethostname())
311
312# distance to nearest wall
313 ywall_s=0.5*(y2d[0:-1,0]+y2d[1:,0])
314 dist_s=yp2d-ywall_s[:,None]
315 ywall_n=0.5*(y2d[0:-1,-1]+y2d[1:,-1])
316 dist_n=ywall_n[:,None] -yp2d
317 dist=np.minimum(dist_s,dist_n)
318
319# west face coordinate
320 xw=0.5*(x2d[0:-1,0:-1]+x2d[0:-1,1:])
321 yw=0.5*(y2d[0:-1,0:-1]+y2d[0:-1,1:])
322
323 del1x=((xw-xp2d)**2+(yw-yp2d)**2)**0.5
324 del2x=((xw-np.roll(xp2d,1,axis=0))**2+(yw-np.roll(yp2d,1,axis=0))**2)**0.5
325 fx=del2x/(del1x+del2x)
326 if cyclic_x:
327 fx[0,:]=0.5
328
329
330# south face coordinate
331 xs=0.5*(x2d[0:-1,0:-1]+x2d[1:,0:-1])
332 ys=0.5*(y2d[0:-1,0:-1]+y2d[1:,0:-1])
333
334 del1y=((xs-xp2d)**2+(ys-yp2d)**2)**0.5
335 del2y=((xs-np.roll(xp2d,1,axis=1))**2+(ys-np.roll(yp2d,1,axis=1))**2)**0.5
336 fy=del2y/(del1y+del2y)
337
338 areawy=np.diff(x2d,axis=1)
339 areawx=-np.diff(y2d,axis=1)
340
341 areasy=-np.diff(x2d,axis=0)
342 areasx=np.diff(y2d,axis=0)
343
344 areaw=(areawx**2+areawy**2)**0.5
345 areas=(areasx**2+areasy**2)**0.5
346
347# volume approaximated as the vector product of two triangles for cells
348 ax=np.diff(x2d,axis=1)
349 ay=np.diff(y2d,axis=1)
350 bx=np.diff(x2d,axis=0)
351 by=np.diff(y2d,axis=0)
352
353 areaz_1=0.5*np.absolute(ax[0:-1,:]*by[:,0:-1]-ay[0:-1,:]*bx[:,0:-1])
354
355 ax=np.diff(x2d,axis=1)
356 ay=np.diff(y2d,axis=1)
357# areaz_2=0.5*np.absolute(ax[1:,:]*by[:,0:-1]-ay[1:,:]*bx[:,0:-1])
358 areaz_2=0.5*np.absolute(ax[1:,:]*by[:,1:]-ay[1:,:]*bx[:,1:])
359
360
361 vol=areaz_1+areaz_2
362
363# coeff at south wall (without viscosity)
364 as_bound=areas[:,0]**2/(0.5*vol[:,0])
365
366# coeff at north wall (without viscosity)
367 an_bound=areas[:,-1]**2/(0.5*vol[:,-1])
368
369# coeff at west wall (without viscosity)
370 aw_bound=areaw[0,:]**2/(0.5*vol[0,:])
371
372 ae_bound=areaw[-1,:]**2/(0.5*vol[-1,:])
373
374 return areaw,areawx,areawy,areas,areasx,areasy,vol,fx,fy,aw_bound,ae_bound,as_bound,an_bound,dist
375
377
378 print('////////////////// Start of input data ////////////////// \n\n\n')
379
380 print('\n\n########### section 1 choice of differencing scheme ###########')
381 print(f"{'scheme: ':<29} {scheme}")
382 print(f"{'scheme_turb: ':<29} {scheme_turb}")
383
384 print('\n\n########### section 2 turbulence models ###########')
385
386 print(f"{'cmu: ':<29} {cmu}")
387 print(f"{'kom: ':<29} {kom}")
388 print(f"{'EARSM: ':<29} {earsm}")
389 print(f"{'PINN: ':<29} {pinn}")
390 if kom:
391 print(f"{'c_omega_1: ':<29} {c_omega_1:.3f}")
392 print(f"{'c_omega_2: ':<29} {c_omega_2}")
393 print(f"{'prand_k: ':<29} {prand_k}")
394 print(f"{'prand_omega: ':<29} {prand_omega}")
395
396 print('\n\n########### section 3 restart/save ###########')
397 print(f"{'restart: ':<29} {restart}")
398 print(f"{'save: ':<29} {save}")
399
400 print('\n\n########### section 4 fluid properties ###########')
401 print(f"{'viscos: ':<29} {viscos:.2e}")
402
403 print('\n\n########### section 5 relaxation factors ###########')
404 print(f"{'urfvis: ':<29} {urfvis}")
405
406 print('\n\n########### section 6 number of iteration and convergence criterira ###########')
407 print(f"{'sormax: ':<29} {sormax}")
408 print(f"{'maxit: ':<29} {maxit}")
409 print(f"{'solver_pp: ':<29} {solver_pp}")
410 print(f"{'solver_vel: ':<29} {solver_vel}")
411 print(f"{'solver_turb: ':<29} {solver_turb}")
412 print(f"{'nsweep_vel: ':<29} {nsweep_vel}")
413 print(f"{'nsweep_pp: ':<29} {nsweep_pp}")
414 print(f"{'nsweep_kom: ':<29} {nsweep_kom}")
415 print(f"{'convergence_limit_u: ':<29} {convergence_limit_u}")
416 print(f"{'convergence_limit_v: ':<29} {convergence_limit_v}")
417 print(f"{'convergence_limit_pp: ':<29} {convergence_limit_pp}")
418 print(f"{'convergence_limit_k: ':<29} {convergence_limit_k}")
419 print(f"{'convergence_limit_om: ':<29} {convergence_limit_om}")
420
421 print('\n\n########### section 7 all variables are printed during the iteration at node ###########')
422 print(f"{'imon: ':<29} {imon}")
423 print(f"{'jmon: ':<29} {jmon}")
424
425
426 print('\n\n########### section 8 time-averaging ###########')
427
428
429 print('\n\n########### section 9 residual scaling parameters ###########')
430 print(f"{'resnorm_p: ':<29} {resnorm_p:.1f}")
431 print(f"{'resnorm_vel: ':<29} {resnorm_vel:.1f}")
432
433
434 print('\n\n########### Section 10 grid and boundary conditions ###########')
435 print(f"{'ni: ':<29} {ni}")
436 print(f"{'nj: ':<29} {nj}")
437 print(f"{'cyclic_x: ':<29} {cyclic_x}")
438 print('\n')
439 print('\n')
440
441 print('------boundary conditions for u')
442 print(f"{' ':<5}{'u_bc_west_type: ':<29} {u_bc_west_type}")
443 print(f"{' ':<5}{'u_bc_east_type: ':<29} {u_bc_east_type}")
444 if u_bc_west_type == 'd':
445 print(f"{' ':<5}{'u_bc_west[0]: ':<29} {u_bc_west[0]}")
446 if u_bc_east_type == 'd':
447 print(f"{' ':<5}{'u_bc_east[0]: ':<29} {u_bc_east[0]}")
448
449
450 print(f"{' ':<5}{'u_bc_south_type: ':<29} {u_bc_south_type}")
451 print(f"{' ':<5}{'u_bc_north_type: ':<29} {u_bc_north_type}")
452
453 if u_bc_south_type == 'd':
454 print(f"{' ':<5}{'u_bc_south[0]: ':<29} {u_bc_south[0]}")
455 if u_bc_north_type == 'd':
456 print(f"{' ':<5}{'u_bc_north[0]: ':<29} {u_bc_north[0]}")
457
458 print('------boundary conditions for v')
459 print(f"{' ':<5}{'v_bc_west_type: ':<29} {v_bc_west_type}")
460 print(f"{' ':<5}{'v_bc_east_type: ':<29} {v_bc_east_type}")
461 if v_bc_west_type == 'd':
462 print(f"{' ':<5}{'v_bc_west[0]: ':<29} {v_bc_west[0]}")
463 if v_bc_east_type == 'd':
464 print(f"{' ':<5}{'v_bc_east[0]: ':<29} {v_bc_east[0]}")
465
466
467 print(f"{' ':<5}{'v_bc_south_type: ':<29} {v_bc_south_type}")
468 print(f"{' ':<5}{'v_bc_north_type: ':<29} {v_bc_north_type}")
469
470 if v_bc_south_type == 'd':
471 print(f"{' ':<5}{'v_bc_south[0]: ':<29} {v_bc_south[0]}")
472 if v_bc_north_type == 'd':
473 print(f"{' ':<5}{'v_bc_north[0]: ':<29} {v_bc_north[0]}")
474
475 print('------boundary conditions for k')
476 print(f"{' ':<5}{'k_bc_west_type: ':<29} {k_bc_west_type}")
477 print(f"{' ':<5}{'k_bc_east_type: ':<29} {k_bc_east_type}")
478 if k_bc_west_type == 'd':
479 print(f"{' ':<5}{'k_bc_west[0]: ':<29} {k_bc_west[0]}")
480 if k_bc_east_type == 'd':
481 print(f"{' ':<5}{'k_bc_east[0]: ':<29} {k_bc_east[0]}")
482
483
484 print(f"{' ':<5}{'k_bc_south_type: ':<29} {k_bc_south_type}")
485 print(f"{' ':<5}{'k_bc_north_type: ':<29} {k_bc_north_type}")
486
487 if k_bc_south_type == 'd':
488 print(f"{' ':<5}{'k_bc_south[0]: ':<29} {k_bc_south[0]}")
489 if k_bc_north_type == 'd':
490 print(f"{' ':<5}{'k_bc_north[0]: ':<29} {k_bc_north[0]}")
491
492 print('------boundary conditions for omega')
493 print(f"{' ':<5}{'om_bc_west_type: ':<29} {om_bc_west_type}")
494 print(f"{' ':<5}{'om_bc_east_type: ':<29} {om_bc_east_type}")
495 if om_bc_west_type == 'd':
496 print(f"{' ':<5}{'om_bc_west[0]: ':<29} {om_bc_west[0]:.1f}")
497 if om_bc_east_type == 'd':
498 print(f"{' ':<5}{'om_bc_east[0]: ':<29} {om_bc_east[0]:.1f}")
499
500 print(f"{' ':<5}{'om_bc_south_type: ':<29} {om_bc_south_type}")
501 print(f"{' ':<5}{'om_bc_north_type: ':<29} {om_bc_north_type}")
502
503 if om_bc_south_type == 'd':
504 print(f"{' ':<5}{'om_bc_south[0]: ':<29} {om_bc_south[0]:.1f}")
505 if om_bc_north_type == 'd':
506 print(f"{' ':<5}{'om_bc_north[0]: ':<29} {om_bc_north[0]:.1f}")
507
508 print('\n\n\n ////////////////// End of input data //////////////////\n\n\n')
509
510 return
511
512def compute_face_phi(phi2d,phi_bc_west,phi_bc_east,phi_bc_south,phi_bc_north,\
513 phi_bc_west_type,phi_bc_east_type,phi_bc_south_type,phi_bc_north_type):
514 import numpy as np
515
516 phi2d_face_w=np.empty((ni+1,nj))
517 phi2d_face_s=np.empty((ni,nj+1))
518 phi2d_face_w[0:-1,:]=fx*phi2d+(1-fx)*np.roll(phi2d,1,axis=0)
519 phi2d_face_s[:,0:-1]=fy*phi2d+(1-fy)*np.roll(phi2d,1,axis=1)
520
521# west boundary
522 phi2d_face_w[0,:]=phi_bc_west
523 if phi_bc_west_type == 'n':
524# neumann
525 phi2d_face_w[0,:]=phi2d[0,:]
526 if cyclic_x:
527 phi2d_face_w[0,:]=0.5*(phi2d[0,:]+phi2d[-1,:])
528
529# east boundary
530 phi2d_face_w[-1,:]=phi_bc_east
531 if phi_bc_east_type == 'n':
532# neumann
533 phi2d_face_w[-1,:]=phi2d[-1,:]
534 phi2d_face_w[-1,:]=phi2d_face_w[-2,:]
535 if cyclic_x:
536 phi2d_face_w[-1,:]=0.5*(phi2d[0,:]+phi2d[-1,:])
537
538
539# south boundary
540 phi2d_face_s[:,0]=phi_bc_south
541 if phi_bc_south_type == 'n':
542# neumann
543 phi2d_face_s[:,0]=phi2d[:,0]
544
545# north boundary
546 phi2d_face_s[:,-1]=phi_bc_north
547 if phi_bc_north_type == 'n':
548# neumann
549 phi2d_face_s[:,-1]=phi2d[:,-1]
550
551 return phi2d_face_w,phi2d_face_s
552
553def dphidx(phi_face_w,phi_face_s):
554
555 phi_w=phi_face_w[0:-1,:]*areawx[0:-1,:]
556 phi_e=-phi_face_w[1:,:]*areawx[1:,:]
557 phi_s=phi_face_s[:,0:-1]*areasx[:,0:-1]
558 phi_n=-phi_face_s[:,1:]*areasx[:,1:]
559 return (phi_w+phi_e+phi_s+phi_n)/vol
560
561def dphidy(phi_face_w,phi_face_s):
562
563 phi_w=phi_face_w[0:-1,:]*areawy[0:-1,:]
564 phi_e=-phi_face_w[1:,:]*areawy[1:,:]
565 phi_s=phi_face_s[:,0:-1]*areasy[:,0:-1]
566 phi_n=-phi_face_s[:,1:]*areasy[:,1:]
567 return (phi_w+phi_e+phi_s+phi_n)/vol
568
569def coeff(convw,convs,vis2d,prand,scheme_local):
570
571 visw=np.zeros((ni+1,nj))
572 viss=np.zeros((ni,nj+1))
573 vis_turb=(vis2d-viscos)/prand
574
575 visw[0:-1,:]=fx*vis_turb+(1-fx)*np.roll(vis_turb,1,axis=0)+viscos
576 viss[:,0:-1]=fy*vis_turb+(1-fy)*np.roll(vis_turb,1,axis=1)+viscos
577
578 volw=np.ones((ni+1,nj))*1e-10
579 vols=np.ones((ni,nj+1))*1e-10
580 volw[1:,:]=0.5*np.roll(vol,-1,axis=0)+0.5*vol
581 diffw=visw[0:-1,:]*areaw[0:-1,:]**2/volw[0:-1,:]
582 vols[:,1:]=0.5*np.roll(vol,-1,axis=1)+0.5*vol
583 diffs=viss[:,0:-1]*areas[:,0:-1]**2/vols[:,0:-1]
584
585 if cyclic_x:
586 visw[0,:]=0.5*(vis_turb[0,:]+vis_turb[-1,:])+viscos
587 diffw[0,:]=visw[0,:]*areaw[0,:]**2/(0.5*(vol[0,:]+vol[-1,:]))
588
589
590 if scheme_local == 'h':
591 if iter == 0:
592 print('hybrid scheme, prand=',prand)
593
594 aw2d=np.maximum(convw[0:-1,:],diffw+(1-fx)*convw[0:-1,:])
595 aw2d=np.maximum(aw2d,0.)
596
597 ae2d=np.maximum(-convw[1:,:],np.roll(diffw,-1,axis=0)-np.roll(fx,-1,axis=0)*convw[1:,:])
598 ae2d=np.maximum(ae2d,0.)
599
600 as2d=np.maximum(convs[:,0:-1],diffs+(1-fy)*convs[:,0:-1])
601 as2d=np.maximum(as2d,0.)
602
603 an2d=np.maximum(-convs[:,1:],np.roll(diffs,-1,axis=1)-np.roll(fy,-1,axis=1)*convs[:,1:])
604 an2d=np.maximum(an2d,0.)
605
606 if scheme_local == 'c':
607 if iter == 0:
608 print('CDS scheme, prand=',prand)
609 aw2d=diffw+(1-fx)*convw[0:-1,:]
610 ae2d=np.roll(diffw,-1,axis=0)-np.roll(fx,-1,axis=0)*convw[1:,:]
611
612 as2d=diffs+(1-fy)*convs[:,0:-1]
613 an2d=np.roll(diffs,-1,axis=1)-np.roll(fy,-1,axis=1)*convs[:,1:]
614
615 if not cyclic_x:
616 aw2d[0,:]=0
617 ae2d[-1,:]=0
618 as2d[:,0]=0
619 an2d[:,-1]=0
620
621 return aw2d,ae2d,as2d,an2d,su2d,sp2d
622
623def bc(su2d,sp2d,phi_bc_west,phi_bc_east,phi_bc_south,phi_bc_north\
624 ,phi_bc_west_type,phi_bc_east_type,phi_bc_south_type,phi_bc_north_type):
625
626 su2d=np.zeros((ni,nj))
627 sp2d=np.zeros((ni,nj))
628
629#south
630 if phi_bc_south_type == 'd':
631 sp2d[:,0]=sp2d[:,0]-viscos*as_bound
632 su2d[:,0]=su2d[:,0]+viscos*as_bound*phi_bc_south
633
634#north
635 if phi_bc_north_type == 'd':
636 sp2d[:,-1]=sp2d[:,-1]-viscos*an_bound
637 su2d[:,-1]=su2d[:,-1]+viscos*an_bound*phi_bc_north
638
639#west
640 if phi_bc_west_type == 'd' and not cyclic_x:
641 sp2d[0,:]=sp2d[0,:]-viscos*aw_bound
642 su2d[0,:]=su2d[0,:]+viscos*aw_bound*phi_bc_west
643#east
644 if phi_bc_east_type == 'd' and not cyclic_x:
645 sp2d[-1,:]=sp2d[-1,:]-viscos*ae_bound
646 su2d[-1,:]=su2d[-1,:]+viscos*ae_bound*phi_bc_east
647
648 return su2d,sp2d
649
650def conv(u2d,v2d,p2d_face_w,p2d_face_s):
651#compute convection
652 u2d_face_w,u2d_face_s=compute_face_phi(u2d,u_bc_west,u_bc_east,u_bc_south,u_bc_north,\
653 u_bc_west_type,u_bc_east_type,u_bc_south_type,u_bc_north_type)
654 v2d_face_w,v2d_face_s=compute_face_phi(v2d,v_bc_west,v_bc_east,v_bc_south,v_bc_north,\
655 v_bc_west_type,v_bc_east_type,v_bc_south_type,v_bc_north_type)
656
657 apw=np.zeros((ni+1,nj))
658 aps=np.zeros((ni,nj+1))
659
660 convw=-u2d_face_w*areawx-v2d_face_w*areawy
661 convs=-u2d_face_s*areasx-v2d_face_s*areasy
662
663#\\\\\\\\\\\\\\\\\ west face
664
665# create ghost cells at east & west boundaries with Neumann b.c.
666 p2d_e=p2d
667 p2d_w=p2d
668# duplicate last row and put it at the end
669 p2d_e=np.insert(p2d_e,-1,p2d_e[-1,:],axis=0)
670# duplicate row 0 and put it before row 0 (west boundary)
671 p2d_w=np.insert(p2d_w,0,p2d_w[0,:],axis=0)
672
673 dp=np.roll(p2d_e,-1,axis=0)-3*p2d_e+3*p2d_w-np.roll(p2d_w,1,axis=0)
674
675# apw[1:,:]=fx*np.roll(ap2d_vel,-1,axis=0)+(1-fx)*ap2d_vel
676 apw[0:-1,:]=fx*ap2d_vel+(1-fx)*np.roll(ap2d_vel,1,axis=0)
677 apw[-1,:]=1e-20
678
679 dvelw=dp*areaw/4/apw
680
681# boundaries (no corrections)
682 dvelw[0,:]=0
683 dvelw[-1,:]=0
684
685 convw=convw+areaw*dvelw
686
687#\\\\\\\\\\\\\\\\\ south face
688# create ghost cells at north & south boundaries with Neumann b.c.
689 p2d_n=p2d
690 p2d_s=p2d
691# duplicate last column and put it at the end
692 p2d_n=np.insert(p2d_n,-1,p2d_n[:,-1],axis=1)
693# duplicate first column and put it before column 0 (south boundary)
694 p2d_s=np.insert(p2d_s,0,p2d_s[:,0],axis=1)
695
696 dp=np.roll(p2d_n,-1,axis=1)-3*p2d_n+3*p2d_s-np.roll(p2d_s,1,axis=1)
697
698# aps[:,1:]=fy*np.roll(ap2d_vel,-1,axis=1)+(1-fy)*ap2d_vel
699 aps[:,0:-1]=fy*ap2d_vel+(1-fy)*np.roll(ap2d_vel,1,axis=1)
700 aps[:,-1]=1e-20
701
702 dvels=dp*areas/4/aps
703
704# boundaries (no corrections)
705 dvels[:,0]=0
706 dvels[:,-1]=0
707
708 convs=convs+areas*dvels
709
710# boundaries
711# west
712 if u_bc_west_type == 'd':
713 convw[0,:]=-u_bc_west*areawx[0,:]-v_bc_west*areawy[0,:]
714# east
715 if u_bc_east_type == 'd':
716 convw[-1,:]=-u_bc_east*areawx[-1,:]-v_bc_east*areawy[-1,:]
717# south
718 if v_bc_south_type == 'd':
719 convs[:,0]=-u_bc_south*areasx[:,0]-v_bc_south*areasy[:,0]
720# north
721 if v_bc_north_type == 'd':
722 convs[:,-1]=-u_bc_north*areasx[:,-1]-v_bc_north*areasy[:,-1]
723
724 return convw,convs
725
726def solve_2d(phi2d,aw2d,ae2d,as2d,an2d,su2d,ap2d,tol_conv,nmax,solver_local):
727 if iter == 0:
728 print('solve_2d called')
729 print('nmax',nmax)
730
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)
736
737 m=ni*nj
738
739 if cyclic_x:
740# A = sparse.diags([ap, -ah[:-1], -al[1:], -an[0:-nk], -as1[nk:], -ae, -aw[nj*nk:],-aw,-ae[nj*nk*(ni-1):]], \
741# [0, 1, -1, nk,-nk, nk*nj, -nk*nj, nj*nk*(ni-1), -nj*nk*(ni-1)], format='csr')
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')
744 else:
745 A = sparse.diags([ap, -an[0:-1], -as1[1:], -ae, -aw[nj:]], [0, 1, -1, nj, -nj], format='csr')
746
747 su=np.matrix.flatten(su2d)
748 phi=np.matrix.flatten(phi2d)
749
750 res_su=np.linalg.norm(su)
751 resid_init=np.linalg.norm(A*phi - su)
752
753 phi_org=phi
754
755# bicg (BIConjugate Gradient)
756# bicgstab (BIConjugate Gradient STABilized)
757# cg (Conjugate Gradient) - symmetric positive definite matrices only
758# cgs (Conjugate Gradient Squared)
759# gmres (Generalized Minimal RESidual)
760# minres (MINimum RESidual)
761# qmr (Quasi
762 abs_tol=1e-10
763 if tol_conv < 0:
764# use absolute convergence criterium
765 abs_tol =abs(tol_conv)*resid_init
766 tol_conv=0
767
768 if solver_local == 'direct':
769 if iter == 0:
770 print('solver in solve_2d: direct solver')
771 info=0
772 resid=np.linalg.norm(A*phi - su)
773 phi = linalg.spsolve(A,su)
774 if solver_local == 'pyamg':
775 if iter == 0:
776 print('solver in solve_2d: pyamg solver')
777 App = pyamg.ruge_stuben_solver(A) # construct the multigrid hierarchy
778 res_amg = []
779 phi = App.solve(su, tol=tol_conv, x0=phi, residuals=res_amg)
780 info=0
781 print('Residual history in pyAMG', ["%0.4e" % i for i in res_amg])
782 if solver_local == 'cgs':
783 if iter == 0:
784 print('solver in solve_2d: cgs')
785 phi,info=linalg.cgs(A,su,x0=phi, rtol=tol_conv, atol=abs_tol, maxiter=nmax) # good
786 if solver_local == 'cg':
787 if iter == 0:
788 print('solver in solve_2d: cg')
789 phi,info=linalg.cg(A,su,x0=phi, rtol=tol_conv, atol=abs_tol, maxiter=nmax) # good
790 if solver_local == 'gmres':
791 if iter == 0:
792 print('solver in solve_2d: gmres')
793 phi,info=linalg.gmres(A,su,x0=phi, rtol=tol_conv, atol=abs_tol, maxiter=nmax) # good
794 if solver_local == 'qmr':
795 if iter == 0:
796 print('solver in solve_2d: qmr')
797 phi,info=linalg.qmr(A,su,x0=phi, rtol=tol_conv, atol=abs_tol, maxiter=nmax) # good
798 if solver_local == 'lgmres':
799 if iter == 0:
800 print('solver in solve_2d: lgmres')
801 phi,info=linalg.lgmres(A,su,x0=phi, rtol=tol_conv, atol=abs_tol, maxiter=nmax) # good
802 if info > 0:
803 print('warning in module solve_2d: convergence in sparse matrix solver not reached')
804
805# compute residual without normalizing with |b|=|su2d|
806 if solver_local != 'direct':
807 resid=np.linalg.norm(A*phi - su)
808
809 delta_phi=np.max(np.abs(phi-phi_org))
810
811 phi2d=np.reshape(phi,(ni,nj))
812 phi2d_org=np.reshape(phi_org,(ni,nj))
813
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}")
817
818# we return the initial residual; otherwise the solution is always satisfied (but the non-linearity is not accounted for)
819 return phi2d,resid_init
820
821def calcu(su2d,sp2d,p2d_face_w,p2d_face_s,uu2d_face_w,uu2d_face_s,uv2d_face_w,uv2d_face_s):
822 if iter == 0:
823 print('calcu called')
824# b.c., sources, coefficients
825
826# presssure gradient
827 dpdx=dphidx(p2d_face_w,p2d_face_s)
828 su2d=su2d-dpdx*vol
829
830# Reynolds stresses (exkl. the part related to vis2d_earsm)
831 if earsm:
832 duudx=dphidx(uu2d_face_w,uu2d_face_s)
833 duvdy=dphidy(uv2d_face_w,uv2d_face_s)
834# phiw=fx(i-1,j,k)*phi(i,j,k,n)+(1.-fx(i-1,j,k))*phi(i-1,j,k,n)
835# phie=fx(i,j,k)*phi(i+1,j,k,n)+(1.-fx(i,j,k))*phi(i,j,k,n)
836# di=(areanx(i,j,k)**2+areany(i,j,k)**2)**0.5
837# dphidi_RANS=(phie-phiw)/di
838 duu_di = (uu2d_face_w[1:,:]-uu2d_face_w[0:-1,:])/areas[:,1:]**0.5
839# phis=fy(i,j-1,k)*phi(i,j,k,n)+(1.-fy(i,j-1,k))*phi(i,j-1,k,n)
840# phin=fy(i,j,k)*phi(i,j+1,k,n)+(1.-fy(i,j,k))*phi(i,j,k,n)
841# dj=(areaex(i,j,k)**2+areaey(i,j,k)**2)**0.5
842# dphidj_RANS=(phin-phis)/dj
843 duv_dj = (uv2d_face_s[:,1:]-uv2d_face_s[:,0:-1])/areaw[1:,:]**0.5
844
845 su2d = su2d-(duudx+duvdy)*vol
846# su2d = su2d-(duu_di+duv_dj)*vol
847
848# modify su & sp
849 su2d,sp2d=modify_u(su2d,sp2d)
850
851 ap2d=aw2d+ae2d+as2d+an2d-sp2d
852
853# under-relaxation
854 ap2d=ap2d/urf_vel
855 su2d=su2d+(1-urf_vel)*ap2d*u2d
856
857 return su2d,sp2d,ap2d
858
859def calcv(su2d,sp2d,p2d_face_w,p2d_face_s,uv2d_face_w,uv2d_face_s,vv2d_face_w,vv2d_face_s):
860 if iter == 0:
861 print('calcv called')
862# b.c., sources, coefficients
863
864# presssure gradient
865 dpdy=dphidy(p2d_face_w,p2d_face_s)
866 su2d=su2d-dpdy*vol
867
868# Reynolds stresses (exkl. the part related to vis2d_earsm)
869 if earsm:
870 duvdx=dphidx(uv2d_face_w,uv2d_face_s)
871 dvvdy=dphidy(vv2d_face_w,vv2d_face_s)
872 duv_di = (uv2d_face_w[1:,:]-uv2d_face_w[0:-1,:])/areas[:,1:]**0.5
873# phis=fy(i,j-1,k)*phi(i,j,k,n)+(1.-fy(i,j-1,k))*phi(i,j-1,k,n)
874# phin=fy(i,j,k)*phi(i,j+1,k,n)+(1.-fy(i,j,k))*phi(i,j,k,n)
875# dj=(areaex(i,j,k)**2+areaey(i,j,k)**2)**0.5
876# dphidj_RANS=(phin-phis)/dj
877 dvv_dj = (vv2d_face_s[:,1:]-vv2d_face_s[:,0:-1])/areaw[1:,:]**0.5
878
879 su2d = su2d-(duvdx+dvvdy)*vol
880# su2d = su2d-(duv_di+dvv_dj)*vol
881
882# modify su & sp
883 su2d,sp2d=modify_v(su2d,sp2d)
884
885 ap2d=aw2d+ae2d+as2d+an2d-sp2d
886
887# under-relaxation
888 ap2d=ap2d/urf_vel
889 su2d=su2d+(1-urf_vel)*ap2d*v2d
890
891# ap2d will be used in calcp; store it as ap2d_vel
892 ap2d_vel=ap2d
893
894 return su2d,sp2d,ap2d,ap2d_vel
895
896def calck(su2d,sp2d,k2d,om2d,vis2d,u2d_face_w,u2d_face_s,v2d_face_w,v2d_face_s):
897# b.c., sources, coefficients
898 if iter == 0:
899 print('calck_kom called')
900
901# production term
902 dudx=dphidx(u2d_face_w,u2d_face_s)
903 dvdx=dphidx(v2d_face_w,v2d_face_s)
904
905 dudy=dphidy(u2d_face_w,u2d_face_s)
906 dvdy=dphidy(v2d_face_w,v2d_face_s)
907
908 if earsm:
909 vist = vis2d_earsm - viscos
910 uu_tot = uu2d - vist*dudx
911 vv_tot = vv2d- vist*dvdy
912 uv_tot = uv2d - vist*(dudy+dvdx)
913 gen = -uu_tot*dudx-uv_tot*(dudy+dvdx)-vv_tot*dvdy
914 su2d=su2d+gen*vol
915 else:
916 gen= (2.*(dudx**2+dvdy**2)+(dudy+dvdx)**2)
917 vist=np.maximum(vis2d-viscos,1e-10)
918 su2d=su2d+vist*gen*vol
919
920# c_k = 1t except when PINN is used, see folder
921# channel-2000-half-channel-PINN
922 sp2d=sp2d-cmu*c_k*om2d*vol
923
924
925# modify su & sp
926 su2d,sp2d=modify_k(su2d,sp2d)
927
928 ap2d=aw2d+ae2d+as2d+an2d-sp2d
929
930# under-relaxation
931 ap2d=ap2d/urf_k
932 su2d=su2d+(1-urf_k)*ap2d*k2d
933
934 return su2d,sp2d,gen,ap2d
935
936def calcom(su2d,sp2d,om2d,gen):
937 if iter == 0:
938 print('calcom called')
939
940#--------production term
941 if earsm:
942 su2d=su2d+c_omega_1*gen*vol*om2d/k2d
943 else:
944 su2d=su2d+c_omega_1*gen*vol
945
946#--------dissipation term
947 sp2d=sp2d-c_omega_2*om2d*vol
948
949# modify su & sp
950 su2d,sp2d=modify_om(su2d,sp2d)
951
952 ap2d=aw2d+ae2d+as2d+an2d-sp2d
953
954# under-relaxation
955 ap2d=ap2d/urf_vel
956 su2d=su2d+(1-urf_omega)*ap2d*om2d
957
958 return su2d,sp2d,ap2d
959
960def calcp(pp2d,ap2d_vel):
961
962 if iter == 0:
963 print('calcp called')
964# b.c., sources, coefficients and under-relaxation for pp2d
965
966 apw=np.zeros((ni+1,nj))
967 aps=np.zeros((ni,nj+1))
968
969 pp2d=0
970#----------simplec: multiply ap by (1-urf)
971 ap2d_vel=np.maximum(ap2d_vel*(1.-urf_vel),1.e-20)
972
973#\\\\\\\\\\\\\\\\ west face
974# visw[0:-1,:,:]=fx*vis_turb+(1-fx)*np.roll(vis_turb,1,axis=0)+viscos
975# viss[:,0:-1,:]=fy*vis_turb+(1-fy)*np.roll(vis_turb,1,axis=1)+viscos
976
977# apw[1:,:]=fx*np.roll(ap2d_vel,-1,axis=0)+(1-fx)*ap2d_vel
978 apw[0:-1,:]=fx*ap2d_vel+(1-fx)*np.roll(ap2d_vel,1,axis=0)
979 if cyclic_x:
980 apw[0,:]=0.5*(ap2d_vel[0,:]+ap2d_vel[-1,:])
981 apw[-1,:]=apw[0,:]
982 else:
983 apw[0,:]=1e-20
984 dw=areawx**2+areawy**2
985 aw2d=dw[0:-1,:]/apw[0:-1,:]
986 ae2d=np.roll(aw2d,-1,axis=0)
987
988#\\\\\\\\\\\\\\\\ south face
989# aps[:,1:]=fy*np.roll(ap2d_vel,-1,axis=1)+(1-fy)*ap2d_vel
990 aps[:,0:-1]=fy*ap2d_vel+(1-fy)*np.roll(ap2d_vel,1,axis=1)
991 aps[:,0]=1e-20
992 ds=areasx**2+areasy**2
993 as2d=ds[:,0:-1]/aps[:,0:-1]
994 an2d=np.roll(as2d,-1,axis=1)
995
996 as2d[:,0]=0
997 an2d[:,-1]=0
998 if not cyclic_x:
999 aw2d[0,:]=0
1000 ae2d[-1,:]=0
1001
1002
1003 ap2d=aw2d+ae2d+as2d+an2d
1004
1005# continuity error
1006# su2d=convw[0:-1,:]-np.roll(convw[0:-1,:],-1,axis=0)+convs[:,0:-1]-np.roll(convs[:,0:-1],-1,axis=1)
1007 su2d=convw[0:-1,:]-convw[1:,:]+convs[:,0:-1]-convs[:,1:]
1008
1009# set pp2d=0 in [0,0] tp make it non-singular
1010 as2d[0,0]=0
1011 an2d[0,0]=0
1012 aw2d[0,0]=0
1013 ae2d[0,0]=0
1014 ap2d[0,0]=1
1015# su2d[0,0]=0
1016
1017 return aw2d,ae2d,as2d,an2d,su2d,ap2d
1018
1019
1020def correct_u_v_p(u2d,v2d,p2d):
1021 if iter == 0:
1022 print('correct_u_v_p called')
1023
1024# correct convections
1025#\\\\\\\\\\\\\ west face
1026 if cyclic_x:
1027# convw[1:-1,:,:]=convw[1:-1,:,:]+aw3d_p[1:,:,:]*(p3d_w[1:-1,:,:]-p3d[1:,:,:])*dtt
1028# convw[0,:,:]=convw[0,:,:]+aw3d_p[0,:,:]*(p3d[-1,:,:]-p3d[0,:,:])*dtt
1029# convw[-1,:,:]=convw[0,:,:]
1030 convw[1:-1,:]=convw[1:-1,:]+aw2d[0:-1,:]*(pp2d[1:,:]-pp2d[0:-1,:])
1031 convw[0,:]=convw[0,:]+aw2d[0,:]*(pp2d[0,:]-pp2d[-1,:])
1032 convw[-1,:]=convw[0,:]
1033 else:
1034 convw[1:-1,:]=convw[1:-1,:]+aw2d[0:-1,:]*(pp2d[1:,:]-pp2d[0:-1,:])
1035
1036#\\\\\\\\\\\\\ south face
1037 convs[:,1:-1]=convs[:,1:-1]+as2d[:,0:-1]*(pp2d[:,1:]-pp2d[:,0:-1])
1038
1039# correct p
1040 p2d=p2d+urf_p*(pp2d-pp2d[0,0])
1041
1042# compute pressure correecion at faces (N.B. p_bc_west,, ... are not used since we impose Neumann b.c., everywhere)
1043 pp2d_face_w,pp2d_face_s=compute_face_phi(pp2d,p_bc_west,p_bc_east,p_bc_south,p_bc_north,\
1044 'n','n','n','n')
1045
1046 dppdx=dphidx(pp2d_face_w,pp2d_face_s)
1047 u2d=u2d-dppdx*vol/ap2d_vel
1048
1049 dppdy=dphidy(pp2d_face_w,pp2d_face_s)
1050 v2d=v2d-dppdy*vol/ap2d_vel
1051
1052
1053 return convw,convs,p2d,u2d,v2d,su2d
1054
1055
1056def vist_kom(vis2d,k2d,om2d):
1057 if iter == 0:
1058 print('vist_kom called')
1059
1060 visold= vis2d
1061 vis2d= k2d/om2d+viscos
1062
1063# modify viscosity
1064 vis2d=modify_vis(vis2d)
1065
1066# under-relax viscosity
1067 vis2d= urfvis*vis2d+(1.-urfvis)*visold
1068
1069 return vis2d
1070
1072 scalar_names = ['pressure']
1073 scalar_variables = [p2d]
1074 scalar_names.append('turb_kin')
1075 scalar_names.append('omega')
1076 scalar_variables.append(k2d)
1077 scalar_variables.append(om2d)
1078
1079 if save_vtk_movie:
1080 file_name = '%s.%d.vtk' % (vtk_file_name, itstep)
1081 else:
1082 file_name = '%s.vtk' % (vtk_file_name)
1083
1084 nk=1
1085 dz=1
1086 f = open(file_name,'w')
1087 f.write('# vtk DataFile Version 3.0\npyCALC-LES Data\nASCII\nDATASET STRUCTURED_GRID\n')
1088 f.write('DIMENSIONS %d %d %d\nPOINTS %d double\n' % (nk+1,nj+1,ni+1,(ni+1)*(nj+1)*(nk+1)))
1089 for i in range(ni+1):
1090 for j in range(nj+1):
1091 for k in range(nk+1):
1092 f.write('%.5f %.5f %.5f\n' % (x2d[i,j],y2d[i,j],dz*k))
1093 f.write('\nCELL_DATA %d\n' % (ni*nj*nk))
1094
1095 f.write('\nVECTORS velocity double\n')
1096 for i in range(ni):
1097 for j in range(nj):
1098 for k in range(nk):
1099 f.write('%.12e %.12e %.12e\n' % (u2d[i,j,k],v2d[i,j,k],w2d[i,j,k]))
1100
1101 for v in range(len(scalar_names)):
1102 var_name = scalar_names[v]
1103 var = scalar_variables[v]
1104 f.write('\nSCALARS %s double 1\nLOOKUP_TABLE default\n' % (var_name))
1105 for i in range(ni):
1106 for j in range(nj):
1107 for k in range(nk):
1108 f.write('%.10e\n' % (var[i,j,k]))
1109 f.close()
1110
1111 print('Flow state save into VTK format to file %s\n' % (file_name))
1112
1114
1115 print('read_restart_data called')
1116
1117 u2d=np.load('u2d_saved.npy')
1118 v2d=np.load('v2d_saved.npy')
1119 p2d=np.load('p2d_saved.npy')
1120 k2d=np.load('k2d_saved.npy')
1121 om2d=np.load('om2d_saved.npy')
1122 vis2d=np.load('vis2d_saved.npy')
1123 ap2d_vel=np.load('ap2d_vel_saved.npy')
1124
1125 return u2d,v2d,p2d,k2d,om2d,vis2d,ap2d_vel
1126
1127def save_data(u2d,v2d,p2d,k2d,om2d,vis2d,ap2d_vel):
1128
1129 print('save_data called')
1130 np.save('u2d_saved', u2d)
1131 np.save('v2d_saved', v2d)
1132 np.save('p2d_saved', p2d)
1133 np.save('k2d_saved', k2d)
1134 np.save('om2d_saved', om2d)
1135 np.save('vis2d_saved', vis2d)
1136 np.save('ap2d_vel_saved', ap2d_vel)
1137 if earsm:
1138 np.save('uu2d_saved', uu2d)
1139 np.save('vv2d_saved', vv2d)
1140 np.save('ww2d_saved', ww2d)
1141 np.save('uv2d_saved', uv2d)
1142 np.save('vis2d_earsm_saved', vis2d_earsm)
1143
1144 return
1145
1146
1147def main():
1148 datax= np.loadtxt("x2d.dat")
1149 x=datax[0:-1]
1150 ni=int(datax[-1])
1151 datay= np.loadtxt("y2d.dat")
1152 y=datay[0:-1]
1153 nj=int(datay[-1])
1154
1155 x2d=np.zeros((ni+1,nj+1))
1156 y2d=np.zeros((ni+1,nj+1))
1157
1158 x2d=np.reshape(x,(ni+1,nj+1))
1159 y2d=np.reshape(y,(ni+1,nj+1))
1160
1161# compute cell centers
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:])
1164
1165# initialize geometric arrays
1166
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))
1182
1183# default values
1184# boundary conditions for u
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)
1189
1190 u_bc_west_type='d'
1191 u_bc_east_type='n'
1192 u_bc_south_type='d'
1193 u_bc_north_type='d'
1194
1195# boundary conditions for v
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)
1200
1201 v_bc_west_type='d'
1202 v_bc_east_type='n'
1203 v_bc_south_type='d'
1204 v_bc_north_type='d'
1205
1206# boundary conditions for p
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)
1211
1212 p_bc_west_type='n'
1213 p_bc_east_type='n'
1214 p_bc_south_type='n'
1215 p_bc_north_type='n'
1216
1217# boundary conditions for k
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)
1222
1223 k_bc_west_type='n'
1224 k_bc_east_type='n'
1225 k_bc_south_type='d'
1226 k_bc_north_type='d'
1227
1228# boundary conditions for omega
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)
1233
1234 om_bc_west_type='d'
1235 om_bc_east_type='n'
1236 om_bc_south_type='d'
1237 om_bc_north_type='d'
1238
1239# boundary conditions for uu
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)
1244
1245 uu_bc_west_type='n'
1246 uu_bc_east_type='n'
1247 uu_bc_south_type='d'
1248 uu_bc_north_type='d'
1249
1250# boundary conditions for uv
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)
1255
1256 uv_bc_west_type='n'
1257 uv_bc_east_type='n'
1258 uv_bc_south_type='d'
1259 uv_bc_north_type='d'
1260
1261# boundary conditions for vv
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)
1266
1267 vv_bc_west_type='n'
1268 vv_bc_east_type='n'
1269 vv_bc_south_type='d'
1270 vv_bc_north_type='d'
1271
1272
1273 cyclic_x=False
1274
1275 vtk=False
1276
1277 earsm = False
1278
1279 pinn = False
1280
1281 setup_case()
1282
1283 print_indata()
1284
1285 areaw,areawx,areawy,areas,areasx,areasy,vol,fx,fy,aw_bound,ae_bound,as_bound,an_bound,dist=init()
1286
1287
1288# initialization
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
1301
1302 fmu2d=np.ones((ni,nj))
1303 gen=np.ones((ni,nj))
1304
1305 convw=np.ones((ni+1,nj))*1e-20
1306 convs=np.ones((ni,nj+1))*1e-20
1307
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
1324
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
1331
1332# comute Delta_max for LES/DES/PANS models
1333 delta_max=np.maximum(vol/areas[:,1:],vol/areaw[1:,:])
1334
1335# the three arrays below are constant except when PINN is used, see folder
1336# channel-2000-half-channel-PINN
1337
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
1341
1342 iter=0
1343
1344# initialize
1345 if pinn:
1346 prand_k, c_k, c_omega_2 = modify_PINN()
1347
1348# initialize
1349 u2d,v2d,k2d,om2d,vis2d,dist=modify_init(u2d,v2d,k2d,om2d,vis2d)
1350
1351# read data from restart
1352 if restart:
1353 u2d,v2d,p2d,k2d,om2d,vis2d,ap2d_vel= read_restart_data()
1354
1355 k2d=np.maximum(k2d,1e-6)
1356
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)
1363
1364
1365 if not cyclic_x:
1366 u_bc_west,v_bc_west,k_bc_west,om_bc_west,u2d_face_w,convw = modify_inlet()
1367
1368 convw,convs=conv(u2d,v2d,p2d_face_w,p2d_face_s)
1369
1370 iter=0
1371
1372 if kom:
1373 urf_temp=urfvis # no under-relaxation
1374 urfvis=1
1375 vis2d=vist_kom(vis2d,k2d,om2d)
1376 urfvis=urf_temp
1377
1378# find max index
1379#sumax=np.max(su2d.flatten())
1380#print('[i,j,k]', np.where(su2d == np.amax(su2d))
1381
1382 residual_u=0
1383 residual_v=0
1384 residual_p=0
1385 residual_k=0
1386 residual_om=0
1387
1388
1389 for iter in range(0,abs(maxit)):
1390
1391 start_time_iter = time.time()
1392# coefficients for velocities
1393 start_time = time.time()
1394# conpute inlet fluc
1395 if iter == 0:
1396 if not cyclic_x:
1397 u_bc_west,v_bc_west,k_bc_west,om_bc_west,u2d_face_w,convw = modify_inlet()
1398 if earsm:
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)
1402 else:
1403 aw2d,ae2d,as2d,an2d,su2d,sp2d=coeff(convw,convs,vis2d,1,scheme)
1404
1405 if earsm:
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)
1412
1413# u2d
1414# boundary conditions for u2d
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)
1418
1419 if maxit > 0:
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}")
1422
1423 start_time = time.time()
1424# v2d
1425# boundary conditions for v2d
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)
1429 if maxit > 0:
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}")
1432
1433 start_time = time.time()
1434# pp2d
1435 convw,convs=conv(u2d,v2d,p2d_face_w,p2d_face_s)
1436 if not cyclic_x:
1437 convw=modify_outlet(convw)
1438 aw2d,ae2d,as2d,an2d,su2d,ap2d=calcp(pp2d,ap2d_vel)
1439 pp2d=np.zeros((ni,nj))
1440 if maxit > 0:
1441 pp2d,residual_p=solve_2d(pp2d,aw2d,ae2d,as2d,an2d,su2d,ap2d,convergence_limit_pp,nsweep_pp,solver_pp)
1442
1443# correct u, v, w, p
1444 convw,convs,p2d,u2d,v2d,su2d= correct_u_v_p(u2d,v2d,p2d)
1445 convw=modify_outlet(convw)
1446
1447# continuity error
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))
1450
1451 print(f"{'time pp: '}{time.time()-start_time:.2e}")
1452
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)
1459
1460 start_time = time.time()
1461
1462 if kom:
1463
1464 if not earsm:
1465 vis2d=vist_kom(vis2d,k2d,om2d)
1466# coefficients
1467 start_time = time.time()
1468 if earsm:
1469 aw2d,ae2d,as2d,an2d,su2d,sp2d=coeff(convw,convs,vis2d_earsm,prand_k,scheme_turb)
1470 else:
1471 aw2d,ae2d,as2d,an2d,su2d,sp2d=coeff(convw,convs,vis2d,prand_k,scheme_turb)
1472
1473# boundary conditions for k2d
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)
1477
1478 aw3d,ae3d,as3d,an3d,ap3d,su3d,sp3d=fix_k()
1479
1480 if maxit > 0:
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}")
1484
1485 start_time = time.time()
1486# omega
1487# boundary conditions for om2d
1488 if earsm:
1489 aw2d,ae2d,as2d,an2d,su2d,sp2d=coeff(convw,convs,vis2d_earsm,prand_omega,scheme_turb)
1490 else:
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)
1495
1496 aw2d,ae2d,as2d,an2d,ap2d,su2d,sp2d=fix_omega()
1497
1498 if maxit > 0:
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)
1501
1502 print(f"{'time omega: '}{time.time()-start_time:.2e}")
1503
1504# scale residuals
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
1510
1511 resmax=np.max([residual_u ,residual_v,residual_pp])
1512
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")
1516
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")
1520
1521
1522
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())
1527
1528 print(f"\n{'---iter: '}{iter:2d}, {'umax: '}{umax:.2e},{'vismax: '}{vismax:.2e}, {'kmin: '}{kmin:.2e}, {'ommin: '}{ommin:.2e}\n")
1529
1530 print(f"{'time one iteration: '}{time.time()-start_time_iter:.2e}")
1531
1532 if resmax < sormax and iter > 0:
1533
1534 break
1535
1536
1537# save data for restart
1538 if save:
1539 save_data(u2d,v2d,p2d,k2d,om2d,vis2d,ap2d_vel)
1540
1541 if vtk:
1542 itstep=ntstep
1543 save_vtk()
1544
1545 print('program reached normal stop')
1546
1547main()
1548
modify_conv(convw, convs)
calcp(pp2d, ap2d_vel)
calcv(su2d, sp2d, p2d_face_w, p2d_face_s, uv2d_face_w, uv2d_face_s, vv2d_face_w, vv2d_face_s)
modify_init(u2d, v2d, k2d, om2d, vis2d)
coeff(convw, convs, vis2d, prand, scheme_local)
bc(su2d, sp2d, phi_bc_west, phi_bc_east, phi_bc_south, phi_bc_north, phi_bc_west_type, phi_bc_east_type, phi_bc_south_type, phi_bc_north_type)
save_data(u2d, v2d, p2d, k2d, om2d, vis2d, ap2d_vel)
dphidy(phi_face_w, phi_face_s)
dphidx(phi_face_w, phi_face_s)
calcu(su2d, sp2d, p2d_face_w, p2d_face_s, uu2d_face_w, uu2d_face_s, uv2d_face_w, uv2d_face_s)
vist_kom(vis2d, k2d, om2d)
compute_face_phi(phi2d, phi_bc_west, phi_bc_east, phi_bc_south, phi_bc_north, phi_bc_west_type, phi_bc_east_type, phi_bc_south_type, phi_bc_north_type)
solve_2d(phi2d, aw2d, ae2d, as2d, an2d, su2d, ap2d, tol_conv, nmax, solver_local)
conv(u2d, v2d, p2d_face_w, p2d_face_s)
calcom(su2d, sp2d, om2d, gen)
correct_u_v_p(u2d, v2d, p2d)
calck(su2d, sp2d, k2d, om2d, vis2d, u2d_face_w, u2d_face_s, v2d_face_w, v2d_face_s)