KTH framework for Nek5000 toolboxes; testing version  0.0.1
rans_komg.f
Go to the documentation of this file.
1  real function rans_komg_mut(ix,iy,iz,iel)
2  include 'SIZE'
3  include 'TSTEP'
4  include 'RANS_KOMG'
5 
6  save ifldla
7  data ifldla /ldimt1/
8 
9  if(ix*iy*iz*iel.eq.1 .and. ifield.le.ifldla) then
10  if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'updating komg_mut'
11  if(ifrans_komg_stndrd) call rans_komg_stndrd_eddy
12  if(ifrans_komg_lowre) call rans_komg_lowre_eddy
13  if(ifrans_komgsst_stndrd) call rans_komgsst_stndrd_eddy
14  if(ifrans_komgsst_lowre) call rans_komgsst_lowre_eddy
15  if(ifrans_komg_stndrd_noreg)call rans_komg_stndrd_eddy_noreg
16  endif
17 
18  ifldla = ifield
19  rans_komg_mut = mut(ix,iy,iz,iel)
20 
21  return
22  end
23 c-----------------------------------------------------------------------
24  real function rans_komg_mutsk(ix,iy,iz,iel)
25  include 'SIZE'
26  include 'TSTEP'
27  include 'RANS_KOMG'
28 
29  rans_komg_mutsk = mutsk(ix,iy,iz,iel)
30 
31  return
32  end
33 c-----------------------------------------------------------------------
34  real function rans_komg_mutso(ix,iy,iz,iel)
35  include 'SIZE'
36  include 'TSTEP'
37  include 'RANS_KOMG'
38 
39  rans_komg_mutso = mutso(ix,iy,iz,iel)
40 
41  return
42  end
43 c-----------------------------------------------------------------------
44  subroutine rans_komg_getcoeffs(coeffs_in)
45  include 'SIZE'
46  include 'TSTEP'
47  include 'RANS_KOMG'
48 
49  real coeffs_in(1)
50 
51  do i=1,ncoeffs
52  coeffs_in(i) = coeffs(i)
53  enddo
54 
55  return
56  end
57 c-----------------------------------------------------------------------
58  real function rans_komg_ksrc(ix,iy,iz,iel)
59  include 'SIZE'
60  include 'TSTEP'
61  include 'RANS_KOMG'
62 
63  logical ifevalsrc
64  data ifevalsrc /.true./
65  common /komgifsrc/ ifevalsrc
66 
67  if(ix*iy*iz*iel.eq.1 .and. ifevalsrc) then
68  if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'updating komg_src'
69  if(ifrans_komg_stndrd) call rans_komg_stndrd_compute
70  if(ifrans_komg_lowre) call rans_komg_lowre_compute
71  if(ifrans_komgsst_stndrd) call rans_komgsst_stndrd_compute
72  if(ifrans_komgsst_lowre) call rans_komgsst_lowre_compute
73  if(ifrans_komg_stndrd_noreg)call rans_komg_stndrd_compute_noreg
74  ifevalsrc = .false.
75  endif
76 
77  if(ifld_k.gt.ifld_omega) ifevalsrc = .true.
78  rans_komg_ksrc = ksrc(ix,iy,iz,iel)
79 
80  return
81  end
82 c-----------------------------------------------------------------------
83  real function rans_komg_omgsrc(ix,iy,iz,iel)
84  include 'SIZE'
85  include 'TSTEP'
86  include 'RANS_KOMG'
87 
88  logical ifevalsrc
89  data ifevalsrc /.true./
90  common /komgifsrc/ ifevalsrc
91 
92  if(ix*iy*iz*iel.eq.1 .and. ifevalsrc) then
93  if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'updating komg_src'
94  if(ifrans_komg_stndrd) call rans_komg_stndrd_compute
95  if(ifrans_komg_lowre) call rans_komg_lowre_compute
96  if(ifrans_komgsst_stndrd) call rans_komgsst_stndrd_compute
97  if(ifrans_komgsst_lowre) call rans_komgsst_lowre_compute
98  if(ifrans_komg_stndrd_noreg)call rans_komg_stndrd_compute_noreg
99  ifevalsrc = .false.
100  endif
101 
102  if(ifld_omega.gt.ifld_k) ifevalsrc = .true.
103  rans_komg_omgsrc = omgsrc(ix,iy,iz,iel)
104 
105  return
106  end
107 c-----------------------------------------------------------------------
108  subroutine rans_komg_stndrd_compute
109 c
110 c Compute RANS source terms and diffusivities on an
111 c element-by-element basis
112 c
113  include 'SIZE'
114  include 'TOTAL'
115  include 'RANS_KOMG'
116 
117  parameter(lxyz=lx1*ly1*lz1)
118 
119  real k_x(lxyz),k_y(lxyz),k_z(lxyz)
120  $ , o_x(lxyz),o_y(lxyz),o_z(lxyz)
121  $ ,omp_x(lxyz), omp_y(lxyz), omp_z(lxyz)
122 
123  real tempv(lxyz), rhoalpk (lxyz)
124  $ ,omwom(lxyz), rhoalpfr(lxyz)
125 
126  real term1(lxyz), term2 (lxyz)
127  $ ,term3(lxyz), term4 (lxyz)
128  $ ,g (lxyz), div (lxyz)
129 
130  common /storesom/ st_mag2(lx1*ly1*lz1,lelv)
131  $ , om_mag2(lx1*ly1*lz1,lelv)
132  $ , oiojsk(lx1*ly1*lz1,lelv)
133  $ , divq(lx1*ly1*lz1,lelv)
134 
135  integer e
136  real k,mu_t,nu_t,mu,nu, kv_min,omeg_max
137 
138  real mu_omeg(lxyz), mu_omegx(lxyz), mu_omegy(lxyz), mu_omegz(lxyz)
139  real extra_src_omega(lxyz)
140 
141 c Turbulent viscosity constants
142  pr_t = coeffs( 1)
143  sigma_k = coeffs( 2)
144  sigma_omega = coeffs( 3)
145 
146 c Low Reynolds number correction constants
147  alpinf_str = coeffs( 4)
148  r_k = coeffs( 5)
149  beta_0 = coeffs( 6)
150  alp0_str = coeffs( 7)
151 
152 c Dissipation of K constants
153  betainf_str = coeffs( 8)
154  alp_inf = coeffs( 9)
155  r_b = coeffs(10)
156  akk = coeffs(11)
157 
158 c Production of omega constants
159  alpha_0 = coeffs(12)
160  r_w = coeffs(13)
161 
162 c Dissipation of omega constants
163 c beta_0 = defined earlier
164 
165  kv_min = coeffs(14)
166  omeg_max = coeffs(15)
167  tiny = coeffs(16)
168 c================================
169 
170  call comp_stom (st_mag2, om_mag2, oiojsk, divq)
171 
172  nome_neg = 0
173  nkey_neg = 0
174  xome_neg = 0.
175  xkey_neg = 0.
176 
177  do e=1,nelv
178 
179  call copy (g, st_mag2(1,e), lxyz)
180  call copy (div, divq(1,e), lxyz)
181 
182  call gradm11(k_x, k_y, k_z, t(1,1,1,1,ifld_k -1),e)
183  call gradm11(omp_x,omp_y,omp_z,t(1,1,1,1,ifld_omega-1),e)
184 
185 c solve for omega_pert
186 
187 c ---------------------
188 c call check_omwall_behavior
189 c ---------------------
190  do i=1,lxyz
191 
192  rho = param(1) ! vtrans(i,1,1,e,1)
193  mu = param(2) ! vdiff (i,1,1,e,1)
194  nu = mu/rho
195 
196 c limits for k, omega
197 
198 c solve for omega_pert
199  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
200  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
201 
202  iflim_omeg = 0 ! limit omega^{prime} 1 limit omega_total
203 
204  if (iflim_omeg.eq.0) then
205  if(t(i,1,1,e,ifld_omega-1).lt.0.0) then
206 c write(*,*) 'Zero OMEG ', t(i,1,1,e,ifld_omega-1)
207  xome_neg = min(xome_neg,t(i,1,1,e,ifld_omega-1))
208  nome_neg = nome_neg + 1
209 c write(*,*) 'Neg OMEG ', nome_neg, xome_neg
210  t(i,1,1,e,ifld_omega-1) =0.01*abs(t(i,1,1,e,ifld_omega-1))
211  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
212  endif
213 
214  if(t(i,1,1,e,ifld_k-1).lt.0.0) then
215 c write(*,*) 'Zero K ', t(i,1,1,e,ifld_k-1)
216  xkey_neg = min(xkey_neg,t(i,1,1,e,ifld_k-1))
217  nkey_neg = nkey_neg + 1
218 c write(*,*) 'Neg KEY ', nkey_neg, xkey_neg
219  t(i,1,1,e,ifld_k-1) = 0.01*abs(t(i,1,1,e,ifld_k-1))
220  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
221  endif
222 
223  elseif(iflim_omeg.eq.1) then
224 
225  if(omega.lt.0.0) then
226 c write(*,*) 'OMEG tot is neg', omega
227  omega = 0.01*abs(omega)
228  endif
229 
230  if(k.lt.0.0) then
231 c write(*,*) 'K is neg', k
232  k = 0.01*abs(k)
233  endif
234 
235  endif
236 
237  expn = -2.0
238  o_x(i)= omp_x(i)+expn * dfdx_omegb(i,1,1,e) *f_omegb(i,1,1,e)
239  o_y(i)= omp_y(i)+expn * dfdy_omegb(i,1,1,e) *f_omegb(i,1,1,e)
240  if(if3d)
241  $ o_z(i)= omp_z(i)+expn * dfdz_omegb(i,1,1,e) *f_omegb(i,1,1,e)
242  omwom(i)= 1.0/(1.0+t(i,1,1,e,ifld_omega-1)/f_omegb(i,1,1,e))
243 
244 c See equations from eqns_k_omega1.pdf from Eq. (3) onwards
245 c Eq.(1) and (2) in eqns_k_omega1.pdf are the governing equations
246 c no source terms Sk or S_w are added
247 
248  st_magn = sqrt(st_mag2(i,e))
249  om_magn = sqrt(om_mag2(i,e))
250  sum_xx = oiojsk(i,e)
251 
252 ! calculate del k * del omega / omega
253 
254  if(if3d)then
255  xk3= (k_x(i)*o_x(i) + k_y(i)*o_y(i) + k_z(i)*o_z(i))
256  $ / (omega**3+tiny)
257  else
258  xk3= (k_x(i)*o_x(i) + k_y(i)*o_y(i))
259  $ / (omega**3+tiny)
260  endif
261 
262  alp_str = alpinf_str
263  mu_t = rho * alp_str*k/(omega + tiny) ! should multiply these by rho!!!
264  rhoalpk(i)= rho*alp_str*k
265 
266  factr = 1.0
267  sigma_omega1 = 1.0/sigma_omega
268  rhoalpfr(i) = expn * rho * alp_str * factr * omwom(i)
269  $ * sigma_omega1
270 
271 c nu_t is kinematic turbulent viscosity
272 c units of k = m2/s2, units of omega = 1/s, nu units = m2/s
273 c set limit for nu_t
274 c nu_t = max(tiny,nu_t)
275 c if(nu_t.gt.5000*mu)nu_t = 5000*mu
276 
277  betai_str = betainf_str
278 
279  if (xk3.le.0)then
280  f_beta_str = 1.0
281  else
282  f_beta_str = (1.0 + 680.0*xk3*xk3)/(1.0 + 400.0*xk3*xk3)
283  endif
284  y_k = rho * betai_str * f_beta_str * k * omega
285 
286 c betai_str = beta_star in 12.5.15 for incompressible flow
287 
288 
289  extra_prod = 0.
290  if(iflomach) extra_prod = twothird*div(i)
291  g_k0= mu_t*g(i) - ( rho*k + mu_t*div(i) )*extra_prod
292  g_k = min(g_k0, 10.*y_k)
293 
294 c g(i) is S**2 in equation sheet
295 
296  ksrc(i,1,1,e) = g_k - y_k
297 
298 c Compute production of omega
299  alpha = (alp_inf/alp_str)
300 
301 c G_w = alpha*alp_str*rho*G_k/mu_t !g(i)
302  g_w0 = alpha*alp_str*rho*(g(i)-(omega+div(i))*extra_prod)
303 
304 c Compute dissipation of omega
305  beta = beta_0
306 c no compressibility correction M < 0.25
307 
308  x_w = abs((sum_xx)/(betainf_str*omega + tiny)**3)
309  f_b = 1.0
310  if(if3d) f_b = (1.0 + 70.0*x_w)/(1.0 + 80.0*x_w)
311 
312  y_w = rho*beta*f_b * omega * omega
313 
314  g_w = min(g_w0, 10.0*y_w)
315  omgsrc(i,1,1,e) = g_w - y_w
316 
317  mut(i,1,1,e) = mu_t
318  mutsk(i,1,1,e) = mu_t / sigma_k
319  mutso(i,1,1,e) = mu_t / sigma_omega
320  enddo
321 
322 c solve for omega_pert
323 
324 c add mut*delsqf
325  sigma_omega1 = 1.0/sigma_omega
326  call copy (mu_omeg,rhoalpk,lxyz)
327  call cmult (mu_omeg,sigma_omega1,lxyz)
328  call col4 (extra_src_omega,mu_omeg ,omwom
329  $ ,delsqf_omegb(1,1,1,e),lxyz)
330 
331 c add mu*delsqf
332  call copy (tempv,delsqf_omegb(1,1,1,e),lxyz)
333  call cmult (tempv,mu,lxyz)
334  call col2 (tempv,f_omegb(1,1,1,e),lxyz)
335  call add2 (extra_src_omega, tempv,lxyz)
336 
337 c Form (1/sigma_w) del_mut * del_omw
338 c form 1: (del_yw/yw) del_k
339  call col3 (term1,dfdx_omegb(1,1,1,e),k_x,lxyz)
340  call addcol3(term1,dfdy_omegb(1,1,1,e),k_y,lxyz)
341  if(if3d)
342  $ call addcol3(term1,dfdz_omegb(1,1,1,e),k_z,lxyz)
343 
344 c form 2: 2(omw/om) (del_omw/yw)^2
345  expm = -expn
346  call col3 (term2,omwom,delfsq_omegb(1,1,1,e), lxyz)
347  call col2 (term2 ,t(1,1,1,e,ifld_k-1),lxyz)
348  call cmult (term2,expm,lxyz)
349  call add3 (tempv, term1, term2, lxyz)
350 
351 c form 3: -(omw/om) k (del_yw/yw) \del_omp/omw
352  call col3 (term3 ,omwom,t(1,1,1,e,ifld_k-1),lxyz)
353  call invcol2(term3 ,f_omegb(1,1,1,e) ,lxyz)
354  call col3 (term4, dfdx_omegb(1,1,1,e),omp_x, lxyz)
355  call addcol3(term4, dfdy_omegb(1,1,1,e),omp_y, lxyz)
356  if(if3d)
357  $ call addcol3(term4, dfdz_omegb(1,1,1,e),omp_z, lxyz)
358  call subcol3(tempv, term3, term4, lxyz)
359 
360  call addcol3(extra_src_omega, rhoalpfr, tempv, lxyz)
361 
362 c add rho v * del_omw
363  call col3 (tempv,dfdx_omegb(1,1,1,e),vx(1,1,1,e),lxyz)
364  call addcol3(tempv,dfdy_omegb(1,1,1,e),vy(1,1,1,e),lxyz)
365  if(if3d)
366  $ call addcol3(tempv,dfdz_omegb(1,1,1,e),vz(1,1,1,e),lxyz)
367  call col2 (tempv,vtrans(1,1,1,e,1),lxyz)
368  call admcol3(extra_src_omega,tempv,f_omegb(1,1,1,e),expm,lxyz)
369 
370  call add2 (omgsrc(1,1,1,e), extra_src_omega ,lxyz)
371 
372  enddo
373 
374  if(loglevel.gt.2) then
375  nome_neg =iglsum(nome_neg,1)
376  nkey_neg =iglsum(nkey_neg,1)
377  xome_neg = glmin(xome_neg,1)
378  xkey_neg = glmin(xkey_neg,1)
379 
380  if(nid.eq.0 .and. nome_neg.gt.0)
381  $ write(*,*) 'Neg Omega ', nome_neg, xome_neg
382  if(nid.eq.0 .and. nkey_neg.gt.0)
383  $ write(*,*) 'Neg TKE ', nkey_neg, xkey_neg
384  endif
385 
386  return
387  end
388 c-----------------------------------------------------------------------
389  subroutine rans_komg_lowre_compute
390 c
391 c Compute RANS source terms and diffusivities on an
392 c element-by-element basis
393 c
394  include 'SIZE'
395  include 'TOTAL'
396  include 'RANS_KOMG'
397 
398  parameter(lxyz=lx1*ly1*lz1)
399 
400  real k_x(lxyz),k_y(lxyz),k_z(lxyz)
401  $ , o_x(lxyz),o_y(lxyz),o_z(lxyz)
402  $ ,omp_x(lxyz), omp_y(lxyz), omp_z(lxyz)
403 
404  real tempv(lxyz), rhoalpk (lxyz)
405  $ ,omwom(lxyz), rhoalpfr(lxyz)
406 
407  real term1(lxyz), term2 (lxyz)
408  $ ,term3(lxyz), term4 (lxyz)
409  $ ,g (lxyz), div (lxyz)
410 
411  common /storesom/ st_mag2(lx1*ly1*lz1,lelv)
412  $ , om_mag2(lx1*ly1*lz1,lelv)
413  $ , oiojsk(lx1*ly1*lz1,lelv)
414  $ , divq(lx1*ly1*lz1,lelv)
415 
416  integer e
417  real k,mu_t,nu_t,mu,nu, kv_min,omeg_max
418 
419  real mu_omeg(lxyz), mu_omegx(lxyz), mu_omegy(lxyz), mu_omegz(lxyz)
420  real extra_src_omega(lxyz)
421 
422 c Turbulent viscosity constants
423  pr_t = coeffs( 1)
424  sigma_k = coeffs( 2)
425  sigma_omega = coeffs( 3)
426 
427 c Low Reynolds number correction constants
428  alpinf_str = coeffs( 4)
429  r_k = coeffs( 5)
430  beta_0 = coeffs( 6)
431  alp0_str = coeffs( 7)
432 
433 c Dissipation of K constants
434  betainf_str = coeffs( 8)
435  alp_inf = coeffs( 9)
436  r_b = coeffs(10)
437  akk = coeffs(11)
438 
439 c Production of omega constants
440  alpha_0 = coeffs(12)
441  r_w = coeffs(13)
442 
443 c Dissipation of omega constants
444 c beta_0 = defined earlier
445 
446  kv_min = coeffs(14)
447  omeg_max = coeffs(15)
448  tiny = coeffs(16)
449 c================================
450 
451  call comp_stom (st_mag2, om_mag2, oiojsk, divq)
452 
453  nome_neg = 0
454  nkey_neg = 0
455  xome_neg = 0.
456  xkey_neg = 0.
457 
458  do e=1,nelv
459 
460  call copy (g, st_mag2(1,e), lxyz)
461  call copy (div, divq(1,e), lxyz)
462 
463  call gradm11(k_x, k_y, k_z, t(1,1,1,1,ifld_k -1),e)
464  call gradm11(omp_x,omp_y,omp_z,t(1,1,1,1,ifld_omega-1),e)
465 
466 c solve for omega_pert
467 
468 c ---------------------
469 c call check_omwall_behavior
470 c ---------------------
471  do i=1,lxyz
472 
473  rho = param(1) ! vtrans(i,1,1,e,1)
474  mu = param(2) ! vdiff (i,1,1,e,1)
475  nu = mu/rho
476 
477 c limits for k, omega
478 
479 c solve for omega_pert
480  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
481  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
482 
483  iflim_omeg = 0 ! limit omega^{prime} 1 limit omega_total
484 
485  if (iflim_omeg.eq.0) then
486  if(t(i,1,1,e,ifld_omega-1).lt.0.0) then
487 c write(*,*) 'Zero OMEG ', t(i,1,1,e,ifld_omega-1)
488  xome_neg = min(xome_neg,t(i,1,1,e,ifld_omega-1))
489  nome_neg = nome_neg + 1
490 c write(*,*) 'Neg OMEG ', nome_neg, xome_neg
491  t(i,1,1,e,ifld_omega-1) =0.01*abs(t(i,1,1,e,ifld_omega-1))
492  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
493  endif
494 
495  if(t(i,1,1,e,ifld_k-1).lt.0.0) then
496 c write(*,*) 'Zero K ', t(i,1,1,e,ifld_k-1)
497  xkey_neg = min(xkey_neg,t(i,1,1,e,ifld_k-1))
498  nkey_neg = nkey_neg + 1
499 c write(*,*) 'Neg KEY ', nkey_neg, xkey_neg
500  t(i,1,1,e,ifld_k-1) = 0.01*abs(t(i,1,1,e,ifld_k-1))
501  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
502  endif
503 
504  elseif(iflim_omeg.eq.1) then
505 
506  if(omega.lt.0.0) then
507  write(*,*) 'OMEG tot is neg', omega
508  omega = 0.01*abs(omega)
509  endif
510 
511  if(k.lt.0.0) then
512  write(*,*) 'K is neg', k
513  k = 0.01*abs(k)
514  endif
515 
516  endif
517 
518  expn = -2.0
519  o_x(i)= omp_x(i)+expn * dfdx_omegb(i,1,1,e) *f_omegb(i,1,1,e)
520  o_y(i)= omp_y(i)+expn * dfdy_omegb(i,1,1,e) *f_omegb(i,1,1,e)
521  if(if3d)
522  $ o_z(i)= omp_z(i)+expn * dfdz_omegb(i,1,1,e) *f_omegb(i,1,1,e)
523  omwom(i)= 1.0/(1.0+t(i,1,1,e,ifld_omega-1)/f_omegb(i,1,1,e))
524 
525 c See equations from eqns_k_omega1.pdf from Eq. (3) onwards
526 c Eq.(1) and (2) in eqns_k_omega1.pdf are the governing equations
527 c no source terms Sk or S_w are added
528 
529  st_magn = sqrt(st_mag2(i,e))
530  om_magn = sqrt(om_mag2(i,e))
531  sum_xx = oiojsk(i,e)
532 
533 ! calculate del k * del omega / omega
534 
535  if(if3d)then
536  xk3= (k_x(i)*o_x(i) + k_y(i)*o_y(i) + k_z(i)*o_z(i))
537  $ / (omega**3+tiny)
538  else
539  xk3= (k_x(i)*o_x(i) + k_y(i)*o_y(i))
540  $ / (omega**3+tiny)
541  endif
542 
543  re_t = rho * k /(mu * omega + tiny)
544  alp_str = alpinf_str * (alp0_str + (re_t/r_k))
545  $ / (1.+(re_t/r_k))
546  mu_t = rho * alp_str*k/(omega + tiny) ! should multiply these by rho!!!
547  rhoalpk(i)= rho*alp_str*k
548 
549  rtld = re_t / r_k
550  funcr = (1.0 - alp0_str) / (alp0_str + rtld)/(1.0 + rtld)
551  factr = (1.0 + rtld * funcr)
552  sigma_omega1 = 1.0/sigma_omega
553  rhoalpfr(i) = expn * rho * alp_str * factr * omwom(i)
554  $ * sigma_omega1
555 
556 c nu_t is kinematic turbulent viscosity
557 c units of k = m2/s2, units of omega = 1/s, nu units = m2/s
558 c set limit for nu_t
559 c nu_t = max(tiny,nu_t)
560 c if(nu_t.gt.5000*mu)nu_t = 5000*mu
561 
562  betai_str = betainf_str * (akk + (re_t/r_b)**4)
563  $ / (1.0 + (re_t/r_b)**4)
564 
565 
566  if (xk3.le.0)then
567  f_beta_str = 1.0
568  else
569  f_beta_str = (1.0 + 680.0*xk3*xk3)/(1.0 + 400.0*xk3*xk3)
570  endif
571  y_k = rho * betai_str * f_beta_str * k * omega
572 
573 c betai_str = beta_star in 12.5.15 for incompressible flow
574 
575  extra_prod = 0.
576  if(iflomach) extra_prod = twothird*div(i)
577  g_k0= mu_t*g(i) - ( rho*k + mu_t*div(i) )*extra_prod
578  g_k = min(g_k0, 10.*y_k)
579 
580 c g(i) is S**2 in equation sheet
581 
582  ksrc(i,1,1,e) = g_k - y_k
583 
584 c Compute production of omega
585  alpha = (alp_inf/alp_str) *
586  $ ((alpha_0 + (re_t/r_w))/(1.0 + (re_t/r_w)))
587 
588 c G_w = alpha*alp_str*rho*G_k/mu_t !g(i)
589  g_w0 = alpha*alp_str*rho*(g(i)-(omega+div(i))*extra_prod)
590 
591 c Compute dissipation of omega
592  beta = beta_0
593 c no compressibility correction M < 0.25
594 
595  x_w = abs((sum_xx)/(betainf_str*omega + tiny)**3)
596  f_b = 1.0
597  if(if3d) f_b = (1.0 + 70.0*x_w)/(1.0 + 80.0*x_w)
598 
599  y_w = rho*beta*f_b * omega * omega
600 
601  g_w = min(g_w0, 10.0*y_w)
602  omgsrc(i,1,1,e) = g_w - y_w
603 
604  mut(i,1,1,e) = mu_t
605  mutsk(i,1,1,e) = mu_t / sigma_k
606  mutso(i,1,1,e) = mu_t / sigma_omega
607  enddo
608 
609 c solve for omega_pert
610 
611 c add mut*delsqf
612  sigma_omega1 = 1.0/sigma_omega
613  call copy (mu_omeg,rhoalpk,lxyz)
614  call cmult (mu_omeg,sigma_omega1,lxyz)
615  call col4 (extra_src_omega,mu_omeg ,omwom
616  $ ,delsqf_omegb(1,1,1,e),lxyz)
617 
618 c add mu*delsqf
619  call copy (tempv,delsqf_omegb(1,1,1,e),lxyz)
620  call cmult (tempv,mu,lxyz)
621  call col2 (tempv,f_omegb(1,1,1,e),lxyz)
622  call add2 (extra_src_omega, tempv,lxyz)
623 
624 c Form (1/sigma_w) del_mut * del_omw
625 c form 1: (del_yw/yw) del_k
626  call col3 (term1,dfdx_omegb(1,1,1,e),k_x,lxyz)
627  call addcol3(term1,dfdy_omegb(1,1,1,e),k_y,lxyz)
628  if(if3d)
629  $ call addcol3(term1,dfdz_omegb(1,1,1,e),k_z,lxyz)
630 
631 c form 2: 2(omw/om) (del_omw/yw)^2
632  expm = -expn
633  call col3 (term2,omwom,delfsq_omegb(1,1,1,e), lxyz)
634  call col2 (term2 ,t(1,1,1,e,ifld_k-1),lxyz)
635  call cmult (term2,expm,lxyz)
636  call add3 (tempv, term1, term2, lxyz)
637 
638 c form 3: -(omw/om) k (del_yw/yw) \del_omp/omw
639  call col3 (term3 ,omwom,t(1,1,1,e,ifld_k-1),lxyz)
640  call invcol2(term3 ,f_omegb(1,1,1,e) ,lxyz)
641  call col3 (term4, dfdx_omegb(1,1,1,e),omp_x, lxyz)
642  call addcol3(term4, dfdy_omegb(1,1,1,e),omp_y, lxyz)
643  if(if3d)
644  $ call addcol3(term4, dfdz_omegb(1,1,1,e),omp_z, lxyz)
645  call subcol3(tempv, term3, term4, lxyz)
646 
647  call addcol3(extra_src_omega, rhoalpfr, tempv, lxyz)
648 
649 c add rho v * del_omw
650  call col3 (tempv,dfdx_omegb(1,1,1,e),vx(1,1,1,e),lxyz)
651  call addcol3(tempv,dfdy_omegb(1,1,1,e),vy(1,1,1,e),lxyz)
652  if(if3d)
653  $ call addcol3(tempv,dfdz_omegb(1,1,1,e),vz(1,1,1,e),lxyz)
654  call col2 (tempv,vtrans(1,1,1,e,1),lxyz)
655  call admcol3(extra_src_omega,tempv,f_omegb(1,1,1,e),expm,lxyz)
656 
657  call add2 (omgsrc(1,1,1,e), extra_src_omega ,lxyz)
658 
659  enddo
660 
661  if(loglevel.gt.2) then
662  nome_neg =iglsum(nome_neg,1)
663  nkey_neg =iglsum(nkey_neg,1)
664  xome_neg = glmin(xome_neg,1)
665  xkey_neg = glmin(xkey_neg,1)
666 
667  if(nid.eq.0 .and. nome_neg.gt.0)
668  $ write(*,*) 'Neg Omega ', nome_neg, xome_neg
669  if(nid.eq.0 .and. nkey_neg.gt.0)
670  $ write(*,*) 'Neg TKE ', nkey_neg, xkey_neg
671  endif
672 
673  return
674  end
675 c-----------------------------------------------------------------------
676  subroutine rans_komgsst_stndrd_compute
677 c
678 c Compute RANS source terms and diffusivities on an
679 c element-by-element basis
680 c
681  include 'SIZE'
682  include 'TOTAL'
683  include 'RANS_KOMG'
684 
685  parameter(lxyz=lx1*ly1*lz1)
686 
687  real k_x(lxyz),k_y(lxyz),k_z(lxyz)
688  $ , o_x(lxyz),o_y(lxyz),o_z(lxyz)
689  $ ,omp_x(lxyz), omp_y(lxyz), omp_z(lxyz)
690 
691  real tempv(lxyz), rhoalpk (lxyz)
692  $ ,omwom(lxyz), rhoalpfr(lxyz)
693 
694  real term1(lxyz), term2 (lxyz)
695  $ ,term3(lxyz), term4 (lxyz)
696  $ ,g (lxyz), div (lxyz)
697 
698  common /storesom/ st_mag2(lx1*ly1*lz1,lelv)
699  $ , om_mag2(lx1*ly1*lz1,lelv)
700  $ , oiojsk(lx1*ly1*lz1,lelv)
701  $ , divq(lx1*ly1*lz1,lelv)
702 
703  integer e
704  real k,mu_t,nu_t,mu,nu, kv_min,omeg_max
705 
706  real mu_omeg(lxyz), mu_omegx(lxyz), mu_omegy(lxyz), mu_omegz(lxyz)
707  real extra_src_omega(lxyz)
708 
709 c====turbulence constants==========
710 
711 c Turbulent viscosity constants
712  pr_t = coeffs( 1)
713  sigk1 = coeffs( 2)
714  sigom1 = coeffs( 3)
715 
716 c Low Reynolds number correction constants
717  alpinf_str = coeffs( 4)
718  r_k = coeffs( 5)
719  beta1 = coeffs( 6)
720  alp0_str = coeffs( 7)
721 
722 c Dissipation of K constants
723  beta_str = coeffs( 8)
724  gamma1 = coeffs( 9)
725  r_b = coeffs(10)
726  akk = coeffs(11)
727 
728 c Production of omega constants
729  alpha_0 = coeffs(12)
730  r_w = coeffs(13)
731 
732 c Dissipation of omega constants
733 c beta_0 = defined earlier
734 
735  kv_min = coeffs(14)
736  omeg_max = coeffs(15)
737  tiny = coeffs(16)
738 
739 c additional SST and k and epsilon constants
740  alp1 = coeffs(17)
741  beta2 = coeffs(18)
742  sigk2 = coeffs(19)
743  sigom2 = coeffs(20)
744  gamma2 = coeffs(21)
745 c================================
746 
747  call comp_stom (st_mag2, om_mag2, oiojsk, divq)
748 
749  nome_neg = 0
750  nkey_neg = 0
751  xome_neg = 0.
752  xkey_neg = 0.
753 
754  do e=1,nelv
755 
756  call copy (g, st_mag2(1,e), lxyz)
757  call copy (div, divq(1,e), lxyz)
758 
759  call gradm11(k_x, k_y, k_z, t(1,1,1,1,ifld_k -1),e)
760  call gradm11(omp_x,omp_y,omp_z,t(1,1,1,1,ifld_omega-1),e)
761 
762 c solve for omega_pert
763 
764 c ---------------------
765 c call check_omwall_behavior
766 c ---------------------
767  do i=1,lxyz
768 
769  rho = param(1) ! vtrans(i,1,1,e,1)
770  mu = param(2) ! vdiff (i,1,1,e,1)
771  nu = mu/rho
772 
773 c limits for k, omega
774 
775 c solve for omega_pert
776  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
777  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
778 
779  iflim_omeg = 0 ! limit omega^{prime} 1 limit omega_total
780 
781  if (iflim_omeg.eq.0) then
782  if(t(i,1,1,e,ifld_omega-1).lt.0.0) then
783 c write(*,*) 'Zero OMEG ', t(i,1,1,e,ifld_omega-1)
784  xome_neg = min(xome_neg,t(i,1,1,e,ifld_omega-1))
785  nome_neg = nome_neg + 1
786 c write(*,*) 'Neg OMEG ', nome_neg, xome_neg
787  t(i,1,1,e,ifld_omega-1) =0.01*abs(t(i,1,1,e,ifld_omega-1))
788  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
789  endif
790 
791  if(t(i,1,1,e,ifld_k-1).lt.0.0) then
792 c write(*,*) 'Zero K ', t(i,1,1,e,ifld_k-1)
793  xkey_neg = min(xkey_neg,t(i,1,1,e,ifld_k-1))
794  nkey_neg = nkey_neg + 1
795 c write(*,*) 'Neg KEY ', nkey_neg, xkey_neg
796  t(i,1,1,e,ifld_k-1) = 0.01*abs(t(i,1,1,e,ifld_k-1))
797  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
798  endif
799 
800  elseif(iflim_omeg.eq.1) then
801 
802  if(omega.lt.0.0) then
803  write(*,*) 'OMEG tot is neg', omega
804  omega = 0.01*abs(omega)
805  endif
806 
807  if(k.lt.0.0) then
808  write(*,*) 'K is neg', k
809  k = 0.01*abs(k)
810  endif
811 
812  endif
813 
814  expn = -2.0
815  o_x(i)= omp_x(i)+expn * dfdx_omegb(i,1,1,e) *f_omegb(i,1,1,e)
816  o_y(i)= omp_y(i)+expn * dfdy_omegb(i,1,1,e) *f_omegb(i,1,1,e)
817  if(if3d)
818  $ o_z(i)= omp_z(i)+expn * dfdz_omegb(i,1,1,e) *f_omegb(i,1,1,e)
819  omwom(i)= 1.0/(1.0+t(i,1,1,e,ifld_omega-1) /f_omegb(i,1,1,e))
820 
821 c See equations from eqns_k_omega1.pdf from Eq. (3) onwards
822 c Eq.(1) and (2) in eqns_k_omega1.pdf are the governing equations
823 c no source terms Sk or S_w are added
824 c ----------
825 
826  st_magn = sqrt(st_mag2(i,e))
827  om_magn = sqrt(om_mag2(i,e))
828 
829 ! calculate del k * del omega / omega
830 
831  if(if3d)then
832  xk = (k_x(i)*o_x(i) + k_y(i)*o_y(i) + k_z(i)*o_z(i))
833  $ / (omega+tiny)
834  else
835  xk = (k_x(i)*o_x(i) + k_y(i)*o_y(i))
836  $ / (omega+tiny)
837  endif
838 
839  rhoalpk(i)= rho*k
840  rhoalpfr(i) = expn * rho * omwom(i) * sigom1
841 
842 ! calculate F2 based on arg2
843 
844  yw = ywd(i,1,1,e)
845  ywm1 = ywdm1(i,1,1,e)
846  ywm2 = ywm1*ywm1
847  arg2_1 = sqrt(abs(k)) * ywm1 / omega / beta_str
848  arg2_2 = 500.0*nu * ywm2 / omega
849  arg2 = max(2.0*arg2_1, arg2_2)
850  fun2 = tanh(arg2 * arg2)
851 
852 ! calculate F1 based on arg1
853 
854  argcd = 2.0 * rho * sigom2 * xk
855  cdkom = max(argcd, 1.0e-10)
856  arg1_1 = max( arg2_1, arg2_2)
857  arg1_2 = 4.0 * rho * sigom2 * k * ywm2 / cdkom
858  arg1 = min( arg1_1, arg1_2)
859  fun1 = tanh(arg1 * arg1 * arg1 * arg1)
860 
861 ! calculate mu_t
862  argn_1 = alp1*(omega + tiny)
863  argn_2 = fun2*st_magn ! this can also be Om_magn
864 
865  mu_t = 0.0
866  denom = max(argn_1, argn_2)
867  if(yw.ne.0) mu_t = rho * alp1 * k / denom
868 
869 ! Compute Y_k = dissipation of k
870 
871  y_k = rho * beta_str * k * omega
872 
873 ! Compute G_k = production of k and limit it to 10*Y_k (the dissipation of k)
874 
875  extra_prod = 0.
876  if(iflomach) extra_prod = twothird*div(i)
877  g_k0= mu_t*g(i) - ( rho*k + mu_t*div(i) )*extra_prod
878  g_k = min(g_k0, 10.*y_k)
879 
880 ! Compute Source term for k
881 
882  ksrc(i,1,1,e) = g_k - y_k
883 
884 c Compute production of omega
885 
886  beta = fun1 * beta1 + (1.0 - fun1) * beta2
887  gamma = fun1 * gamma1 + (1.0 - fun1) * gamma2
888  sigk = fun1 * sigk1 + (1.0 - fun1) * sigk2
889  sigom = fun1 * sigom1 + (1.0 - fun1) * sigom2
890 
891 ! Compute production of omega
892 
893  g_w0= rho * gamma * (g(i)-(div(i)+denom/alp1)*extra_prod)
894 
895 ! Compute dissipation of omega
896 
897  y_w = rho * beta * omega * omega
898 
899  g_w = min(g_w0, 10.0*y_w)
900 
901 ! Compute additional SST term for omega
902 
903  s_w = (1.0 - fun1) * argcd
904 
905 ! Compute Source term for omega
906 
907  omgsrc(i,1,1,e) = g_w - y_w + s_w
908  mut(i,1,1,e) = mu_t
909  mutsk(i,1,1,e) = mu_t * sigk
910  mutso(i,1,1,e) = mu_t * sigom
911 
912  enddo
913 
914 c solve for omega_pert
915 
916 c add mut*delsqf
917  call copy (mu_omeg,rhoalpk,lxyz)
918  call cmult (mu_omeg,sigom1, lxyz)
919  call col4 (extra_src_omega,mu_omeg ,omwom
920  $ ,delsqf_omegb(1,1,1,e),lxyz)
921 
922 c add mu*delsqf
923  call copy (tempv,delsqf_omegb(1,1,1,e),lxyz)
924  call cmult (tempv,mu,lxyz)
925  call col2 (tempv,f_omegb(1,1,1,e),lxyz)
926  call add2 (extra_src_omega, tempv,lxyz)
927 
928 c Form (1/sigma_w) del_mut * del_omw
929 c form 1: (del_yw/yw) del_k
930  call col3 (term1,dfdx_omegb(1,1,1,e),k_x,lxyz)
931  call addcol3(term1,dfdy_omegb(1,1,1,e),k_y,lxyz)
932  if(if3d)
933  $ call addcol3(term1,dfdz_omegb(1,1,1,e),k_z,lxyz)
934 
935 c form 2: 2(omw/om) (del_omw/yw)^2
936  expm = -expn
937  call col3 (term2,omwom,delfsq_omegb(1,1,1,e), lxyz)
938  call col2 (term2 ,t(1,1,1,e,ifld_k-1),lxyz)
939  call cmult (term2,expm,lxyz)
940  call add3 (tempv, term1, term2, lxyz)
941 
942 c form 3: -(omw/om) k (del_yw/yw) \del_omp/omw
943  call col3 (term3 ,omwom,t(1,1,1,e,ifld_k-1),lxyz)
944  call invcol2(term3 ,f_omegb(1,1,1,e) ,lxyz)
945  call col3 (term4, dfdx_omegb(1,1,1,e),omp_x, lxyz)
946  call addcol3(term4, dfdy_omegb(1,1,1,e),omp_y, lxyz)
947  if(if3d)
948  $ call addcol3(term4, dfdz_omegb(1,1,1,e),omp_z, lxyz)
949  call subcol3(tempv, term3, term4, lxyz)
950 
951  call addcol3(extra_src_omega, rhoalpfr, tempv, lxyz)
952 
953 c add rho v * del_omw
954  call col3 (tempv,dfdx_omegb(1,1,1,e),vx(1,1,1,e),lxyz)
955  call addcol3(tempv,dfdy_omegb(1,1,1,e),vy(1,1,1,e),lxyz)
956  if(if3d)
957  $ call addcol3(tempv,dfdz_omegb(1,1,1,e),vz(1,1,1,e),lxyz)
958  call col2 (tempv,vtrans(1,1,1,e,1),lxyz)
959  call admcol3(extra_src_omega,tempv,f_omegb(1,1,1,e),expm,lxyz)
960 
961  call add2 (omgsrc(1,1,1,e), extra_src_omega ,lxyz)
962 
963  enddo
964 
965  if(loglevel.gt.2) then
966  nome_neg =iglsum(nome_neg,1)
967  nkey_neg =iglsum(nkey_neg,1)
968  xome_neg = glmin(xome_neg,1)
969  xkey_neg = glmin(xkey_neg,1)
970 
971  if(nid.eq.0 .and. nome_neg.gt.0)
972  $ write(*,*) 'Neg Omega ', nome_neg, xome_neg
973  if(nid.eq.0 .and. nkey_neg.gt.0)
974  $ write(*,*) 'Neg TKE ', nkey_neg, xkey_neg
975  endif
976 
977  return
978  end
979 c-----------------------------------------------------------------------
980  subroutine rans_komgsst_lowre_compute
981 c
982 c Compute RANS source terms and diffusivities on an
983 c element-by-element basis
984 c
985  include 'SIZE'
986  include 'TOTAL'
987  include 'RANS_KOMG'
988 
989  parameter(lxyz=lx1*ly1*lz1)
990 
991  real k_x(lxyz),k_y(lxyz),k_z(lxyz)
992  $ , o_x(lxyz),o_y(lxyz),o_z(lxyz)
993  $ ,omp_x(lxyz), omp_y(lxyz), omp_z(lxyz)
994 
995  real tempv(lxyz), rhoalpk (lxyz)
996  $ ,omwom(lxyz), rhoalpfr(lxyz)
997 
998  real term1(lxyz), term2 (lxyz)
999  $ ,term3(lxyz), term4 (lxyz)
1000  $ ,g (lxyz), div (lxyz)
1001 
1002  common /storesom/ st_mag2(lx1*ly1*lz1,lelv)
1003  $ , om_mag2(lx1*ly1*lz1,lelv)
1004  $ , oiojsk(lx1*ly1*lz1,lelv)
1005  $ , divq(lx1*ly1*lz1,lelv)
1006 
1007  integer e
1008  real k,mu_t,nu_t,mu,nu, kv_min,omeg_max
1009 
1010  real mu_omeg(lxyz), mu_omegx(lxyz), mu_omegy(lxyz), mu_omegz(lxyz)
1011  real extra_src_omega(lxyz)
1012 
1013 c====turbulence constants==========
1014 
1015 c Turbulent viscosity constants
1016  pr_t = coeffs( 1)
1017  sigk1 = coeffs( 2)
1018  sigom1 = coeffs( 3)
1019 
1020 c Low Reynolds number correction constants
1021  alpinf_str = coeffs( 4)
1022  r_k = coeffs( 5)
1023  beta1 = coeffs( 6)
1024  alp0_str = coeffs( 7)
1025 
1026 c Dissipation of K constants
1027  beta_str = coeffs( 8)
1028  alp_inf = coeffs( 9)
1029  r_b = coeffs(10)
1030  akk = coeffs(11)
1031 
1032 c Production of omega constants
1033  alpha_0 = coeffs(12)
1034  r_w = coeffs(13)
1035 
1036 c Dissipation of omega constants
1037 c beta_0 = defined earlier
1038 
1039  kv_min = coeffs(14)
1040  omeg_max = coeffs(15)
1041  tiny = coeffs(16)
1042 
1043 c additional SST and k and epsilon constants
1044  alp1 = coeffs(17)
1045  beta2 = coeffs(18)
1046  sigk2 = coeffs(19)
1047  sigom2 = coeffs(20)
1048  gamma2 = coeffs(21)
1049 c================================
1050 
1051  call comp_stom (st_mag2, om_mag2, oiojsk, divq)
1052 
1053  nome_neg = 0
1054  nkey_neg = 0
1055  xome_neg = 0.
1056  xkey_neg = 0.
1057 
1058  do e=1,nelv
1059 
1060  call copy (g, st_mag2(1,e), lxyz)
1061  call copy (div, divq(1,e), lxyz)
1062 
1063  call gradm11(k_x, k_y, k_z, t(1,1,1,1,ifld_k -1),e)
1064  call gradm11(omp_x,omp_y,omp_z,t(1,1,1,1,ifld_omega-1),e)
1065 
1066 c solve for omega_pert
1067 
1068 c ---------------------
1069 c call check_omwall_behavior
1070 c ---------------------
1071  do i=1,lxyz
1072 
1073  rho = param(1) ! vtrans(i,1,1,e,1)
1074  mu = param(2) ! vdiff (i,1,1,e,1)
1075  nu = mu/rho
1076 
1077 c limits for k, omega
1078 
1079 c solve for omega_pert
1080  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
1081  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
1082 
1083  iflim_omeg = 0 ! limit omega^{prime} 1 limit omega_total
1084 
1085  if (iflim_omeg.eq.0) then
1086  if(t(i,1,1,e,ifld_omega-1).lt.0.0) then
1087 c write(*,*) 'Zero OMEG ', t(i,1,1,e,ifld_omega-1)
1088  xome_neg = min(xome_neg,t(i,1,1,e,ifld_omega-1))
1089  nome_neg = nome_neg + 1
1090 c write(*,*) 'Neg OMEG ', nome_neg, xome_neg
1091  t(i,1,1,e,ifld_omega-1) =0.01*abs(t(i,1,1,e,ifld_omega-1))
1092  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
1093  endif
1094 
1095  if(t(i,1,1,e,ifld_k-1).lt.0.0) then
1096 c write(*,*) 'Zero K ', t(i,1,1,e,ifld_k-1)
1097  xkey_neg = min(xkey_neg,t(i,1,1,e,ifld_k-1))
1098  nkey_neg = nkey_neg + 1
1099 c write(*,*) 'Neg KEY ', nkey_neg, xkey_neg
1100  t(i,1,1,e,ifld_k-1) = 0.01*abs(t(i,1,1,e,ifld_k-1))
1101  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
1102  endif
1103 
1104  elseif(iflim_omeg.eq.1) then
1105 
1106  if(omega.lt.0.0) then
1107  write(*,*) 'OMEG tot is neg', omega
1108  omega = 0.01*abs(omega)
1109  endif
1110 
1111  if(k.lt.0.0) then
1112  write(*,*) 'K is neg', k
1113  k = 0.01*abs(k)
1114  endif
1115 
1116  endif
1117 
1118  expn = -2.0
1119  o_x(i)= omp_x(i)+expn * dfdx_omegb(i,1,1,e) *f_omegb(i,1,1,e)
1120  o_y(i)= omp_y(i)+expn * dfdy_omegb(i,1,1,e) *f_omegb(i,1,1,e)
1121  if(if3d)
1122  $ o_z(i)= omp_z(i)+expn * dfdz_omegb(i,1,1,e) *f_omegb(i,1,1,e)
1123  omwom(i)= 1.0/(1.0+t(i,1,1,e,ifld_omega-1)/f_omegb(i,1,1,e))
1124 
1125 c See equations from eqns_k_omega1.pdf from Eq. (3) onwards
1126 c Eq.(1) and (2) in eqns_k_omega1.pdf are the governing equations
1127 c no source terms Sk or S_w are added
1128 c ----------
1129 
1130  st_magn = sqrt(st_mag2(i,e))
1131  om_magn = sqrt(om_mag2(i,e))
1132  sum_xx = oiojsk(i,e)
1133 
1134 ! calculate del k * del omega / omega
1135 
1136  if(if3d)then
1137  xk = (k_x(i)*o_x(i) + k_y(i)*o_y(i) + k_z(i)*o_z(i))
1138  $ / (omega+tiny)
1139  xk3= (k_x(i)*o_x(i) + k_y(i)*o_y(i) + k_z(i)*o_z(i))
1140  $ / (omega**3+tiny)
1141  else
1142  xk = (k_x(i)*o_x(i) + k_y(i)*o_y(i))
1143  $ / (omega+tiny)
1144  xk3= (k_x(i)*o_x(i) + k_y(i)*o_y(i))
1145  $ / (omega**3+tiny)
1146  endif
1147 
1148 c ------------ begin low Re part of k-omega -------------------
1149 
1150  re_t = rho * k /(mu * omega + tiny)
1151  alp_str = alpinf_str * (alp0_str + (re_t/r_k))
1152  $ / (1.+(re_t/r_k))
1153  rhoalpk(i)= rho*alp_str*k
1154 
1155  rtld = re_t / r_k
1156  funcr = (1.0 - alp0_str) / (alp0_str + rtld)/(1.0 + rtld)
1157  factr = (1.0 + rtld * funcr)
1158  rhoalpfr(i) = expn * rho * alp_str * factr * omwom(i) * sigom1
1159 
1160 c nu_t is kinematic turbulent viscosity
1161 c units of k = m2/s2, units of omega = 1/s, nu units = m2/s
1162 c set limit for nu_t
1163 c nu_t = max(tiny,nu_t)
1164 c if(nu_t.gt.5000*mu)nu_t = 5000*mu
1165 
1166  betai_str = beta_str * (akk + (re_t/r_b)**4)
1167  $ / (1.0 + (re_t/r_b)**4)
1168 
1169  if (xk3.le.0)then
1170  f_beta_str = 1.0
1171  else
1172  f_beta_str = (1.0 + 680.0*xk3*xk3)/(1.0 + 400.0*xk3*xk3)
1173  endif
1174 
1175 c Compute production of omega
1176  alpha = (alp_inf/alp_str) *
1177  $ ((alpha_0 + (re_t/r_w))/(1.0 + (re_t/r_w)))
1178 
1179  gammai = alpha*alp_str
1180 
1181 c Compute dissipation of omega
1182 
1183  x_w = abs((sum_xx)/(beta_str*omega + tiny)**3)
1184  f_b = 1.0
1185  if(if3d) f_b = (1.0 + 70.0*x_w)/(1.0 + 80.0*x_w)
1186  betai= beta1 * f_b
1187 
1188 c ------------ end low Re part of k-omega -------------------
1189 
1190 ! calculate F2 based on arg2
1191 
1192  yw = ywd(i,1,1,e)
1193  ywm1 = ywdm1(i,1,1,e)
1194  ywm2 = ywm1*ywm1
1195  arg2_1 = sqrt(abs(k)) * ywm1 / omega / beta_str
1196  arg2_2 = 500.0*nu * ywm2 / omega
1197  arg2 = max(2.0*arg2_1, arg2_2)
1198  fun2 = tanh(arg2 * arg2)
1199 
1200 ! calculate F1 based on arg1
1201 
1202  argcd = 2.0 * rho * sigom2 * xk
1203  cdkom = max(argcd, 1.0e-10)
1204  arg1_1 = max( arg2_1, arg2_2)
1205  arg1_2 = 4.0 * rho * sigom2 * k * ywm2 / cdkom
1206  arg1 = min( arg1_1, arg1_2)
1207  fun1 = tanh(arg1 * arg1 * arg1 * arg1)
1208 
1209 ! calculate mu_t
1210  argn_1 = alp1*(omega + tiny) / alp_str
1211  argn_2 = fun2*st_magn ! this can also be St_magn
1212 
1213  mu_t = 0.0
1214  denom = max(argn_1, argn_2)
1215  if(yw.ne.0) mu_t = rho * alp1 * k / denom
1216 
1217 ! Compute Y_k = dissipation of k
1218 
1219 c Y_k = rho * beta_str * k * omega
1220  y_k = rho * betai_str * f_beta_str * k * omega
1221 
1222 ! Compute G_k = production of k and limit it to 10*Y_k (the dissipation of k)
1223 
1224  extra_prod = 0.
1225  if(iflomach) extra_prod = twothird*div(i)
1226  g_k0= mu_t*g(i) - ( rho*k + mu_t*div(i) )*extra_prod
1227  g_k = min(g_k0, 10.*y_k)
1228 
1229 ! Compute Source term for k
1230 
1231  ksrc(i,1,1,e) = g_k - y_k
1232 
1233 c Compute production of omega
1234 
1235  beta = fun1 * betai + (1.0 - fun1) * beta2
1236  gamma = fun1 * gammai + (1.0 - fun1) * gamma2
1237  sigk = fun1 * sigk1 + (1.0 - fun1) * sigk2
1238  sigom = fun1 * sigom1 + (1.0 - fun1) * sigom2
1239 
1240 ! Compute production of omega
1241 
1242  g_w0= rho * gamma * (g(i)-(div(i)+denom/alp1)*extra_prod)
1243 
1244 ! Compute dissipation of omega
1245 
1246  y_w = rho * beta * omega * omega
1247 
1248  g_w = min(g_w0, 10.0*y_w)
1249 
1250 ! Compute additional SST term for omega
1251 
1252  s_w = (1.0 - fun1) * argcd
1253 
1254 ! Compute Source term for omega
1255 
1256  omgsrc(i,1,1,e) = g_w - y_w + s_w
1257  mut(i,1,1,e) = mu_t
1258  mutsk(i,1,1,e) = mu_t * sigk
1259  mutso(i,1,1,e) = mu_t * sigom
1260 
1261  enddo
1262 
1263 c solve for omega_pert
1264 
1265 c add mut*delsqf
1266  call copy (mu_omeg,rhoalpk,lxyz)
1267  call cmult (mu_omeg,sigom1, lxyz)
1268  call col4 (extra_src_omega,mu_omeg ,omwom
1269  $ ,delsqf_omegb(1,1,1,e),lxyz)
1270 
1271 c add mu*delsqf
1272  call copy (tempv,delsqf_omegb(1,1,1,e),lxyz)
1273  call cmult (tempv,mu,lxyz)
1274  call col2 (tempv,f_omegb(1,1,1,e),lxyz)
1275  call add2 (extra_src_omega, tempv,lxyz)
1276 
1277 c Form (1/sigma_w) del_mut * del_omw
1278 c form 1: (del_yw/yw) del_k
1279  call col3 (term1,dfdx_omegb(1,1,1,e),k_x,lxyz)
1280  call addcol3(term1,dfdy_omegb(1,1,1,e),k_y,lxyz)
1281  if(if3d)
1282  $ call addcol3(term1,dfdz_omegb(1,1,1,e),k_z,lxyz)
1283 
1284 c form 2: 2(omw/om) (del_omw/yw)^2
1285  expm = -expn
1286  call col3 (term2,omwom,delfsq_omegb(1,1,1,e), lxyz)
1287  call col2 (term2 ,t(1,1,1,e,ifld_k-1),lxyz)
1288  call cmult (term2,expm,lxyz)
1289  call add3 (tempv, term1, term2, lxyz)
1290 
1291 c form 3: -(omw/om) k (del_yw/yw) \del_omp/omw
1292  call col3 (term3 ,omwom,t(1,1,1,e,ifld_k-1),lxyz)
1293  call invcol2(term3 ,f_omegb(1,1,1,e) ,lxyz)
1294  call col3 (term4, dfdx_omegb(1,1,1,e),omp_x, lxyz)
1295  call addcol3(term4, dfdy_omegb(1,1,1,e),omp_y, lxyz)
1296  if(if3d)
1297  $ call addcol3(term4, dfdz_omegb(1,1,1,e),omp_z, lxyz)
1298  call subcol3(tempv, term3, term4, lxyz)
1299 
1300  call addcol3(extra_src_omega, rhoalpfr, tempv, lxyz)
1301 
1302 c add rho v * del_omw
1303  call col3 (tempv,dfdx_omegb(1,1,1,e),vx(1,1,1,e),lxyz)
1304  call addcol3(tempv,dfdy_omegb(1,1,1,e),vy(1,1,1,e),lxyz)
1305  if(if3d)
1306  $ call addcol3(tempv,dfdz_omegb(1,1,1,e),vz(1,1,1,e),lxyz)
1307  call col2 (tempv,vtrans(1,1,1,e,1),lxyz)
1308  call admcol3(extra_src_omega,tempv,f_omegb(1,1,1,e),expm,lxyz)
1309 
1310  call add2 (omgsrc(1,1,1,e), extra_src_omega ,lxyz)
1311 
1312  enddo
1313 
1314  if(loglevel.gt.2) then
1315  nome_neg =iglsum(nome_neg,1)
1316  nkey_neg =iglsum(nkey_neg,1)
1317  xome_neg = glmin(xome_neg,1)
1318  xkey_neg = glmin(xkey_neg,1)
1319 
1320  if(nid.eq.0 .and. nome_neg.gt.0)
1321  $ write(*,*) 'Neg Omega ', nome_neg, xome_neg
1322  if(nid.eq.0 .and. nkey_neg.gt.0)
1323  $ write(*,*) 'Neg TKE ', nkey_neg, xkey_neg
1324  endif
1325 
1326  return
1327  end
1328 c-----------------------------------------------------------------------
1329  subroutine rans_komg_stndrd_compute_noreg
1330 c
1331 c Compute RANS source terms and diffusivities on an
1332 c element-by-element basis
1333 c
1334  include 'SIZE'
1335  include 'TOTAL'
1336  include 'RANS_KOMG'
1337 
1338  parameter(lxyz=lx1*ly1*lz1)
1339 
1340  real k_x(lxyz),k_y(lxyz),k_z(lxyz)
1341  $ , o_x(lxyz),o_y(lxyz),o_z(lxyz)
1342  $ ,omp_x(lxyz), omp_y(lxyz), omp_z(lxyz)
1343 
1344  real tempv(lxyz), rhoalpk (lxyz)
1345  $ ,omwom(lxyz), rhoalpfr(lxyz)
1346 
1347  real term1(lxyz), term2 (lxyz)
1348  $ ,term3(lxyz), term4 (lxyz)
1349  $ ,g (lxyz), div (lxyz)
1350 
1351  common /storesom/ st_mag2(lx1*ly1*lz1,lelv)
1352  $ , om_mag2(lx1*ly1*lz1,lelv)
1353  $ , oiojsk(lx1*ly1*lz1,lelv)
1354  $ , divq(lx1*ly1*lz1,lelv)
1355 
1356  integer e
1357  real k,mu_t,nu_t,mu,nu, kv_min,omeg_max
1358 
1359  real mu_omeg(lxyz), mu_omegx(lxyz), mu_omegy(lxyz), mu_omegz(lxyz)
1360  real extra_src_omega(lxyz)
1361 
1362 c Turbulent viscosity constants
1363  pr_t = coeffs( 1)
1364  sigma_k = coeffs( 2)
1365  sigma_omega = coeffs( 3)
1366 
1367 c Low Reynolds number correction constants
1368  alpinf_str = coeffs( 4)
1369  r_k = coeffs( 5)
1370  beta_0 = coeffs( 6)
1371  alp0_str = coeffs( 7)
1372 
1373 c Dissipation of K constants
1374  betainf_str = coeffs( 8)
1375  alp_inf = coeffs( 9)
1376  r_b = coeffs(10)
1377  akk = coeffs(11)
1378 
1379 c Production of omega constants
1380  alpha_0 = coeffs(12)
1381  r_w = coeffs(13)
1382 
1383 c Dissipation of omega constants
1384 c beta_0 = defined earlier
1385 
1386  kv_min = coeffs(14)
1387  omeg_max = coeffs(15)
1388  tiny = coeffs(16)
1389 c================================
1390 
1391  call comp_stom (st_mag2, om_mag2, oiojsk, divq)
1392 
1393  nome_neg = 0
1394  nkey_neg = 0
1395  xome_neg = 0.
1396  xkey_neg = 0.
1397 
1398  do e=1,nelv
1399 
1400  call copy (g, st_mag2(1,e), lxyz)
1401  call copy (div, divq(1,e), lxyz)
1402 c call rzero (div, lxyz)
1403 
1404  call gradm11(k_x, k_y, k_z, t(1,1,1,1,ifld_k -1),e)
1405  call gradm11(omp_x,omp_y,omp_z,t(1,1,1,1,ifld_omega-1),e)
1406 
1407 c solve for omega_pert
1408 
1409 c ---------------------
1410 c call check_omwall_behavior
1411 c ---------------------
1412  do i=1,lxyz
1413 
1414  rho = param(1) ! vtrans(i,1,1,e,1)
1415  mu = param(2) ! vdiff (i,1,1,e,1)
1416  nu = mu/rho
1417 
1418 c limits for k, omega
1419 
1420 c solve for omega_pert
1421  omega = t(i,1,1,e,ifld_omega-1) ! Current k & omega values
1422  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
1423 
1424  iflim_omeg = 0 ! limit omega^{prime} 1 limit omega_total
1425 
1426  if (iflim_omeg.eq.0) then
1427  if(t(i,1,1,e,ifld_omega-1).lt.0.0) then
1428 c write(*,*) 'Zero OMEG ', t(i,1,1,e,ifld_omega-1)
1429  xome_neg = min(xome_neg,t(i,1,1,e,ifld_omega-1))
1430  nome_neg = nome_neg + 1
1431 c write(*,*) 'Neg OMEG ', nome_neg, xome_neg
1432  t(i,1,1,e,ifld_omega-1) =0.01*abs(t(i,1,1,e,ifld_omega-1))
1433  omega = t(i,1,1,e,ifld_omega-1) ! Current k & omega values
1434  endif
1435 
1436  if(t(i,1,1,e,ifld_k-1).lt.0.0) then
1437 c write(*,*) 'Zero K ', t(i,1,1,e,ifld_k-1)
1438  xkey_neg = min(xkey_neg,t(i,1,1,e,ifld_k-1))
1439  nkey_neg = nkey_neg + 1
1440 c write(*,*) 'Neg KEY ', nkey_neg, xkey_neg
1441  t(i,1,1,e,ifld_k-1) = 0.01*abs(t(i,1,1,e,ifld_k-1))
1442  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
1443  endif
1444 
1445  elseif(iflim_omeg.eq.1) then
1446 
1447  if(omega.lt.0.0) then
1448  write(*,*) 'OMEG tot is neg', omega
1449  omega = 0.01*abs(omega)
1450  endif
1451 
1452  if(k.lt.0.0) then
1453  write(*,*) 'K is neg', k
1454  k = 0.01*abs(k)
1455  endif
1456 
1457  endif
1458 
1459  expn = -2.0
1460  o_x(i)= omp_x(i)
1461  o_y(i)= omp_y(i)
1462  if(if3d)
1463  $ o_z(i)= omp_z(i)
1464 
1465 c See equations from eqns_k_omega1.pdf from Eq. (3) onwards
1466 c Eq.(1) and (2) in eqns_k_omega1.pdf are the governing equations
1467 c no source terms Sk or S_w are added
1468 
1469  if(st_mag2(i,e).lt.0.) write(*,*) ' St_mag2', i, e, st_mag2
1470  if(om_mag2(i,e).lt.0.) write(*,*) ' Om_mag2', i, e, om_mag2
1471  st_magn = sqrt(st_mag2(i,e))
1472  om_magn = sqrt(om_mag2(i,e))
1473  sum_xx = oiojsk(i,e)
1474 c if(sum_xx .ne.0.) write(*,*) ' sum_xx ', i, e, sum_xx
1475 
1476 ! calculate del k * del omega / omega
1477 
1478  if(if3d)then
1479  xk3= (k_x(i)*o_x(i) + k_y(i)*o_y(i) + k_z(i)*o_z(i))
1480  $ / (omega**3+tiny)
1481  else
1482  xk3= (k_x(i)*o_x(i) + k_y(i)*o_y(i))
1483  $ / (omega**3+tiny)
1484  endif
1485 
1486  alp_str = alpinf_str
1487  mu_t = rho * alp_str*k/(omega + tiny) ! should multiply these by rho!!!
1488 
1489 c nu_t is kinematic turbulent viscosity
1490 c units of k = m2/s2, units of omega = 1/s, nu units = m2/s
1491 c set limit for nu_t
1492 c nu_t = max(tiny,nu_t)
1493 c if(nu_t.gt.5000*mu)nu_t = 5000*mu
1494 
1495  betai_str = betainf_str
1496 
1497  if (xk3.le.0)then
1498  f_beta_str = 1.0
1499  else
1500  f_beta_str = (1.0 + 680.0*xk3*xk3)/(1.0 + 400.0*xk3*xk3)
1501  endif
1502  y_k = rho * betai_str * f_beta_str * k * omega
1503 
1504 c betai_str = beta_star in 12.5.15 for incompressible flow
1505 
1506  extra_prod = 0.
1507  if(iflomach) extra_prod = twothird*div(i)
1508  g_k0= mu_t*g(i) - ( rho*k + mu_t*div(i) )*extra_prod
1509  g_k = min(g_k0, 10.*y_k)
1510 
1511 c g(i) is S**2 in equation sheet
1512 
1513  ksrc(i,1,1,e) = g_k - y_k
1514 
1515 c Compute production of omega
1516  alpha = (alp_inf/alp_str)
1517 
1518 c G_w0 = alpha*alp_str*rho*g(i)
1519  g_w0 = alpha*alp_str*rho*(g(i)-(omega+div(i))*extra_prod)
1520 
1521 c Compute dissipation of omega
1522  beta = beta_0
1523 c no compressibility correction M < 0.25
1524 
1525  x_w = abs((sum_xx)/(betainf_str*omega + tiny)**3)
1526  f_b = 1.0
1527  if(if3d) f_b = (1.0 + 70.0*x_w)/(1.0 + 80.0*x_w)
1528 
1529  y_w = rho*beta*f_b * omega * omega
1530 
1531  g_w = min(g_w0, 10.0*y_w)
1532  omgsrc(i,1,1,e) = g_w - y_w
1533 
1534  mut(i,1,1,e) = mu_t
1535  mutsk(i,1,1,e) = mu_t / sigma_k
1536  mutso(i,1,1,e) = mu_t / sigma_omega
1537  enddo
1538 
1539  enddo
1540 
1541  if(loglevel.gt.2) then
1542  nome_neg =iglsum(nome_neg,1)
1543  nkey_neg =iglsum(nkey_neg,1)
1544  xome_neg = glmin(xome_neg,1)
1545  xkey_neg = glmin(xkey_neg,1)
1546 
1547  if(nid.eq.0 .and. nome_neg.gt.0)
1548  $ write(*,*) 'Neg Omega ', nome_neg, xome_neg
1549  if(nid.eq.0 .and. nkey_neg.gt.0)
1550  $ write(*,*) 'Neg TKE ', nkey_neg, xkey_neg
1551  endif
1552 
1553  return
1554  end
1555 c-----------------------------------------------------------------------
1556  subroutine rans_komg_init(ifld_k_in,ifld_omega_in,ifcoeffs
1557  $ ,coeffs_in,wall_id,ywd_in,model_id)
1558 c
1559 c Initialize values ifld_omega & ifld_k for RANS k-omega turbulence
1560 c modeling
1561 c
1562  include 'SIZE'
1563  include 'TOTAL'
1564  include 'RANS_KOMG'
1565 
1566  real w1,w2,w3,w4,w5
1567  common /scrns/
1568  & w1(lx1*ly1*lz1*lelv)
1569  &,w2(lx1*ly1*lz1*lelv)
1570  &,w3(lx1*ly1*lz1*lelv)
1571  &,w4(lx1*ly1*lz1*lelv)
1572  &,w5(lx1*ly1*lz1*lelv)
1573 
1574  integer n,wall_id,ifld_mx
1575  real coeffs_in(1),ywd_in(1)
1576  logical ifcoeffs
1577 
1578  character*3 bcw
1579  character*36 mname(6)
1580 
1581  data mname
1582  &/'regularized standard k-omega '
1583  &,'regularized low-Re k-omega '
1584  &,'regularized standard k-omega SST '
1585  &,'regularized low-Re k-omega SST '
1586  &,'non-regularized standard k-omega '
1587  &,'non-regularized standard k-omega SST'/
1588 
1589  n=nx1*ny1*nz1*nelv
1590 
1591  if(nid.eq.0) write(6,*) 'init RANS model'
1592 
1593  if(iflomach) then
1594  if(nid.eq.0) write(6,*)
1595  & "ERROR: K-OMEGA NOT SUPPORTED WITH LOW MACH FORMULATION"
1596  call exitt
1597  endif
1598 
1599  ifrans_komg_stndrd = .false.
1600  ifrans_komg_lowre = .false.
1601  ifrans_komgsst_stndrd = .false.
1602  ifrans_komgsst_lowre = .false.
1603  ifrans_komg_stndrd_noreg = .false.
1604  if(model_id .eq.0) ifrans_komg_stndrd = .true.
1605  if(model_id .eq.1) ifrans_komg_lowre = .true.
1606  if(model_id .eq.2) ifrans_komgsst_stndrd = .true.
1607  if(model_id .eq.3) ifrans_komgsst_lowre = .true.
1608  if(model_id .eq.4) ifrans_komg_stndrd_noreg = .true.
1609 
1610  if(nid.eq.0) write(*,'(a,a)')
1611  & ' model: ',mname(model_id+1)
1612  ifld_k = ifld_k_in
1613  ifld_omega = ifld_omega_in
1614  ifld_mx=max(ifld_k,ifld_omega)
1615  if (ifld_mx.gt.ldimt1)
1616  $ call exitti('nflds gt ldimt+1, recompile with ldimt > ',
1617  $ ifld_mx+1)
1618 
1619 ! specify k-omega model coefficients
1620 
1621 c if(ncoeffs_in.lt.ncoeffs)
1622 c $ call exitti('dim of user provided komg coeffs array
1623 c $ should be >=$',ncoeffs)
1624 
1625  if(ifcoeffs) then
1626  do i=1,ncoeffs
1627  coeffs(i) =coeffs_in(i)
1628  enddo
1629  else
1630  if(ifrans_komg_stndrd .or. ifrans_komg_lowre .or.
1631  $ ifrans_komg_stndrd_noreg)call rans_komg_set_defaultcoeffs
1632  if(ifrans_komgsst_stndrd .or. ifrans_komgsst_lowre)
1633  $ call rans_komgsst_set_defaultcoeffs
1634  endif
1635 
1636 c solve for omega_pert
1637  if(wall_id.eq.0) then
1638  if(nid.eq.0) write(6,*) ' user supplied wall distance'
1639  call copy(ywd,ywd_in,n)
1640  else
1641  bcw = 'W '
1642  ifld = 1
1643  if(nid.eq.0) write(6,*) 'BC for distance ',bcw
1644  if(wall_id.eq.1) call cheap_dist(ywd,ifld,bcw)
1645  if(wall_id.eq.2) call distf(ywd,ifld,bcw,w1,w2,w3,w4,w5)
1646  call copy(ywd_in,ywd,n)
1647  endif
1648 
1649  call rans_komg_omegabase
1650 
1651  if(nid.eq.0) write(6,*) 'done :: init RANS'
1652 
1653  return
1654  end
1655 c-----------------------------------------------------------------------
1656  subroutine rans_komg_omegabase
1657 c
1658 c Compute Omega base solution which diverges at walls
1659 c
1660  include 'SIZE'
1661  include 'TOTAL'
1662  include 'RANS_KOMG'
1663 
1664  parameter(nprof_max = 10000)
1665  integer i,j,k,e
1666  real nu
1667 
1668  real kv_min,omeg_max
1669 
1670  omeg_max = coeffs(15)
1671  beta0 = coeffs(6)
1672  nu = param(2)/param(1)
1673  ntot1 = nx1*ny1*nz1*nelv
1674 
1675  betainf_str = coeffs(11)
1676  yw_min = glmin(ywd,ntot1)
1677 
1678  cfcon = 6.0 * nu / beta0 ! 2.0 * nu0 / betainf_str !
1679  expn = -2.0
1680 
1681  call gradm1 (dfdx_omegb,dfdy_omegb,dfdz_omegb, ywd)
1682  call opcolv (dfdx_omegb,dfdy_omegb,dfdz_omegb, bm1)
1683  call opdssum(dfdx_omegb,dfdy_omegb,dfdz_omegb)
1684  call opcolv (dfdx_omegb,dfdy_omegb,dfdz_omegb,binvm1)
1685 
1686 c call gradm1 (dudx, dudy ,dudz, dfdx_omegb)
1687 c call copy (delsqf_omegb, dudx, ntot1)
1688 c call gradm1 (dudx, dudy ,dudz, dfdy_omegb)
1689 c call add2 (delsqf_omegb, dudy, ntot1)
1690 c call gradm1 (dudx, dudy ,dudz, dfdz_omegb)
1691 c call add2 (delsqf_omegb, dudz, ntot1)
1692 
1693  toll = 1.0e-08
1694 
1695 c call vdot3 (delfsq_omegb,dfdx_omegb,dfdy_omegb,dfdz_omegb
1696 c $ ,dfdx_omegb,dfdy_omegb,dfdz_omegb,ntot1)
1697 
1698  call rzero (delsqf_omegb, ntot1)
1699  call rone (delfsq_omegb, ntot1)
1700 
1701  do e = 1,nelv
1702  do k = 1,nz1
1703  do j = 1,ny1
1704  do i = 1,nx1
1705 
1706  ieg = lglel(e)
1707  yw = ywd(i,j,k,e) ! 1.0 - abs(y)
1708  ywmin=sqrt(cfcon/omeg_max)
1709  if(yw.gt.ywmin) then
1710  ywm1 = 1.0 /yw
1711  else
1712  if(yw.ne.0) write(*,'(a,3G14.7,4I5)')
1713  $ 'ywmin and yw ',ywmin,yw,omeg_max, i, j, k, ieg
1714  ywm1 = 1.0 /ywmin
1715  endif
1716  ywdm1(i,j,k,e) = ywm1
1717  ywm2 = ywm1*ywm1
1718  ywm3 = ywm2*ywm1
1719  ywm4 = ywm2*ywm2
1720 
1721  f_omegb(i,j,k,e) = cfcon * ywm2
1722  delfpart = expn * ywm1
1723  dfdx_omegb(i,j,k,e) = dfdx_omegb(i,j,k,e) * ywm1
1724  dfdy_omegb(i,j,k,e) = dfdy_omegb(i,j,k,e) * ywm1
1725  dfdz_omegb(i,j,k,e) = dfdz_omegb(i,j,k,e) * ywm1
1726  delsqfpart = expn * (expn - 1.0) * ywm2
1727  delsqf_omegb(i,j,k,e) =(delsqfpart * delfsq_omegb(i,j,k,e)
1728  $ + delfpart * delsqf_omegb(i,j,k,e))
1729  delfsq_omegb(i,j,k,e) = delfsq_omegb(i,j,k,e)* ywm2
1730 
1731  enddo
1732  enddo
1733  enddo
1734  enddo
1735 
1736  return
1737  end
1738 c-----------------------------------------------------------------------
1739  subroutine rans_komg_set_defaultcoeffs
1740 c
1741  include 'SIZE'
1742  include 'RANS_KOMG'
1743 
1744 c ====various problem-specific turbulence constants
1745 c omeg_max = value of omega on the walls
1746 c kv_min = value of K on the walls
1747 c Pr_t is the turbulent prandtl number
1748 
1749  vkappa = 0.41
1750 
1751 c Turbulent viscosity constants
1752  pr_t = 0.85
1753  coeffs( 1) = pr_t
1754  sigma_k = 2.0
1755  coeffs( 2) = sigma_k
1756  sigma_omega = 2.0
1757  coeffs( 3) = sigma_omega
1758 
1759 c Low Reynolds number correction constants
1760 
1761 c Production of K constants
1762  alpinf_str = 1.0
1763  coeffs( 4) = alpinf_str
1764  r_k = 6.0
1765  coeffs( 5) = r_k
1766  beta_0 = 0.072 ! should be 0.075 for SST
1767  coeffs( 6) = beta_0
1768  alp0_str = beta_0/3.0
1769  coeffs( 7) = alp0_str
1770 
1771 c Dissipation of K constants
1772  betainf_str = 0.09
1773  coeffs( 8) = betainf_str
1774  alp_inf = 0.52
1775 c alp_inf = beta_0/betainf_str
1776 c $ - vkappa**2/sqrt(betainf_str)/sigma_omega ! should be 0.52 for k-omega
1777  coeffs( 9) = alp_inf
1778  r_b = 8.0
1779  coeffs(10) = r_b
1780  akk = 4.0/15.0
1781  coeffs(11) = akk
1782 
1783 c Production of omega constants
1784  alpha_0 = 1.0/9.0
1785  coeffs(12) = alpha_0
1786  r_w = 2.95
1787  coeffs(13) = r_w
1788 
1789 c Dissipation of omega constants
1790 c beta_0 = defined earlier
1791 
1792  kv_min = 0.0
1793  coeffs(14) = kv_min
1794  omeg_max = 2.0e8 ! 400.0 Lan
1795  coeffs(15) = omeg_max
1796  tiny = 1.0e-20
1797  coeffs(16) = tiny
1798 
1799  return
1800  end
1801 c-----------------------------------------------------------------------
1802  subroutine rans_komgsst_set_defaultcoeffs
1803 c
1804  include 'SIZE'
1805  include 'RANS_KOMG'
1806 
1807 c ====various problem-specific turbulence constants
1808 c omeg_max = value of omega on the walls
1809 c kv_min = value of K on the walls
1810 c Pr_t is the turbulent prandtl number
1811 
1812  vkappa = 0.41
1813 
1814 c Turbulent viscosity constants
1815  pr_t = 0.85
1816  coeffs( 1) = pr_t
1817  sigk1 = 0.5 ! should be 0.85 for SST
1818  coeffs( 2) = sigk1
1819  sigom1 = 0.5
1820  coeffs( 3) = sigom1
1821 
1822 c Low Reynolds number correction constants
1823 
1824 c Production of K constants
1825  alpinf_str = 1.0
1826  coeffs( 4) = alpinf_str
1827  r_k = 6.0
1828  coeffs( 5) = r_k
1829  beta_0 = 0.072 ! should be 0.075 for SST
1830  coeffs( 6) = beta_0
1831  alp0_str = beta_0/3.0
1832  coeffs( 7) = alp0_str
1833 
1834 c Dissipation of K constants
1835  betainf_str = 0.09
1836  coeffs( 8) = betainf_str
1837  alp_inf = beta_0/betainf_str
1838  $ - sigom1 * vkappa**2/sqrt(betainf_str) ! should be 0.52 for k-omega
1839 c alp_inf = 0.52
1840  coeffs( 9) = alp_inf
1841  r_b = 8.0
1842  coeffs(10) = r_b
1843  akk = 4.0/15.0
1844  coeffs(11) = akk
1845 
1846 c Production of omega constants
1847  alpha_0 = 1.0/9.0
1848  coeffs(12) = alpha_0
1849  r_w = 2.95
1850  coeffs(13) = r_w
1851 
1852 c Dissipation of omega constants
1853 c beta_0 = defined earlier
1854 
1855  kv_min = 0.0
1856  coeffs(14) = kv_min
1857  omeg_max = 2.0e8 ! 400.0
1858  coeffs(15) = omeg_max
1859  tiny = 1.0e-20
1860  coeffs(16) = tiny
1861 
1862 c additional SST and k and epsilon constants
1863  alp1 = 0.31
1864  coeffs(17) = alp1
1865  beta2 = 0.0828
1866  coeffs(18) = beta2
1867  sigk2 = 1.0
1868  coeffs(19) = sigk2
1869  sigom2 = 0.856
1870  coeffs(20) = sigom2
1871  gamma2 = beta2/betainf_str
1872  $ - sigom2 * vkappa**2/sqrt(betainf_str) ! should be 0.44 for k-epsilon
1873 c gamma2 = 0.44
1874  coeffs(21) = gamma2
1875 
1876  return
1877  end
1878 c-----------------------------------------------------------------------
1879  subroutine rans_komg_stndrd_eddy
1880 c
1881 c Compute RANS source terms and diffusivities on an
1882 c element-by-element basis
1883 c
1884  include 'SIZE'
1885  include 'TOTAL'
1886  include 'RANS_KOMG'
1887 
1888  parameter(lxyz=lx1*ly1*lz1)
1889 
1890  integer e
1891  real k,mu_t,nu_t,mu,nu, kv_min,omeg_max
1892 
1893 c Turbulent viscosity constants
1894  pr_t = coeffs( 1)
1895  sigma_k = coeffs( 2)
1896  sigma_omega = coeffs( 3)
1897 
1898 c Low Reynolds number correction constants
1899  alpinf_str = coeffs( 4)
1900  r_k = coeffs( 5)
1901  beta_0 = coeffs( 6)
1902  alp0_str = coeffs( 7)
1903 
1904 c Dissipation of K constants
1905  betainf_str = coeffs( 8)
1906  alp_inf = coeffs( 9)
1907  r_b = coeffs(10)
1908  akk = coeffs(11)
1909 
1910 c Production of omega constants
1911  alpha_0 = coeffs(12)
1912  r_w = coeffs(13)
1913 
1914 c Dissipation of omega constants
1915 c beta_0 = defined earlier
1916 
1917  kv_min = coeffs(14)
1918  omeg_max = coeffs(15)
1919  tiny = coeffs(16)
1920 c================================
1921 
1922  nome_neg = 0
1923  nkey_neg = 0
1924  xome_neg = 0.
1925  xkey_neg = 0.
1926 
1927  do e=1,nelv
1928  do i=1,lxyz
1929 
1930  rho = param(1) ! vtrans(i,1,1,e,1)
1931  mu = param(2) ! vdiff (i,1,1,e,1)
1932  nu = mu/rho
1933 
1934 c limits for k, omega
1935 
1936 c solve for omega_pert
1937 
1938  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
1939  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
1940 
1941  iflim_omeg = 0 ! limit omega^{prime} 1 limit omega_total
1942 
1943  if (iflim_omeg.eq.0) then
1944  if(t(i,1,1,e,ifld_omega-1).lt.0.0) then
1945 c write(*,*) 'Zero OMEG ', t(i,1,1,e,ifld_omega-1)
1946  xome_neg = min(xome_neg,t(i,1,1,e,ifld_omega-1))
1947  nome_neg = nome_neg + 1
1948 c write(*,*) 'Neg OMEG ', nome_neg, xome_neg
1949  t(i,1,1,e,ifld_omega-1) =0.01*abs(t(i,1,1,e,ifld_omega-1))
1950  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
1951  endif
1952 
1953  if(t(i,1,1,e,ifld_k-1).lt.0.0) then
1954 c write(*,*) 'Zero K ', t(i,1,1,e,ifld_k-1)
1955  xkey_neg = min(xkey_neg,t(i,1,1,e,ifld_k-1))
1956  nkey_neg = nkey_neg + 1
1957 c write(*,*) 'Neg KEY ', nkey_neg, xkey_neg
1958  t(i,1,1,e,ifld_k-1) = 0.01*abs(t(i,1,1,e,ifld_k-1))
1959  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
1960  endif
1961 
1962  elseif(iflim_omeg.eq.1) then
1963 
1964  if(omega.lt.0.0) then
1965  write(*,*) 'OMEG tot is neg', omega
1966  omega = 0.01*abs(omega)
1967  endif
1968 
1969  if(k.lt.0.0) then
1970  write(*,*) 'K is neg', k
1971  k = 0.01*abs(k)
1972  endif
1973 
1974  endif
1975 
1976  alp_str = alpinf_str
1977  mu_t = rho * alp_str*k/(omega + tiny) ! should multiply these by rho!!!
1978 
1979  mut(i,1,1,e) = mu_t
1980  mutsk(i,1,1,e) = mu_t / sigma_k
1981  mutso(i,1,1,e) = mu_t / sigma_omega
1982  enddo
1983 
1984  enddo
1985 
1986  if(loglevel.gt.2) then
1987  nome_neg =iglsum(nome_neg,1)
1988  nkey_neg =iglsum(nkey_neg,1)
1989  xome_neg = glmin(xome_neg,1)
1990  xkey_neg = glmin(xkey_neg,1)
1991 
1992  if(nid.eq.0 .and. nome_neg.gt.0)
1993  $ write(*,*) 'Neg Omega ', nome_neg, xome_neg
1994  if(nid.eq.0 .and. nkey_neg.gt.0)
1995  $ write(*,*) 'Neg TKE ', nkey_neg, xkey_neg
1996  endif
1997 
1998  return
1999  end
2000 c-----------------------------------------------------------------------
2001  subroutine rans_komg_lowre_eddy
2002 c
2003 c Compute RANS source terms and diffusivities on an
2004 c element-by-element basis
2005 c
2006  include 'SIZE'
2007  include 'TOTAL'
2008  include 'RANS_KOMG'
2009 
2010  parameter(lxyz=lx1*ly1*lz1)
2011 
2012  integer e
2013  real k,mu_t,nu_t,mu,nu, kv_min,omeg_max
2014 
2015 c Turbulent viscosity constants
2016  pr_t = coeffs( 1)
2017  sigma_k = coeffs( 2)
2018  sigma_omega = coeffs( 3)
2019 
2020 c Low Reynolds number correction constants
2021  alpinf_str = coeffs( 4)
2022  r_k = coeffs( 5)
2023  beta_0 = coeffs( 6)
2024  alp0_str = coeffs( 7)
2025 
2026 c Dissipation of K constants
2027  betainf_str = coeffs( 8)
2028  alp_inf = coeffs( 9)
2029  r_b = coeffs(10)
2030  akk = coeffs(11)
2031 
2032 c Production of omega constants
2033  alpha_0 = coeffs(12)
2034  r_w = coeffs(13)
2035 
2036 c Dissipation of omega constants
2037 c beta_0 = defined earlier
2038 
2039  kv_min = coeffs(14)
2040  omeg_max = coeffs(15)
2041  tiny = coeffs(16)
2042 c================================
2043 
2044  nome_neg = 0
2045  nkey_neg = 0
2046  xome_neg = 0.
2047  xkey_neg = 0.
2048 
2049  do e=1,nelv
2050 
2051  do i=1,lxyz
2052 
2053  rho = param(1) ! vtrans(i,1,1,e,1)
2054  mu = param(2) ! vdiff (i,1,1,e,1)
2055  nu = mu/rho
2056 
2057 c limits for k, omega
2058 
2059 c solve for omega_pert
2060  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
2061  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
2062 
2063  iflim_omeg = 0 ! limit omega^{prime} 1 limit omega_total
2064 
2065  if (iflim_omeg.eq.0) then
2066  if(t(i,1,1,e,ifld_omega-1).lt.0.0) then
2067 c write(*,*) 'Zero OMEG ', t(i,1,1,e,ifld_omega-1)
2068  xome_neg = min(xome_neg,t(i,1,1,e,ifld_omega-1))
2069  nome_neg = nome_neg + 1
2070 c write(*,*) 'Neg OMEG ', nome_neg, xome_neg
2071  t(i,1,1,e,ifld_omega-1) =0.01*abs(t(i,1,1,e,ifld_omega-1))
2072  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
2073  endif
2074 
2075  if(t(i,1,1,e,ifld_k-1).lt.0.0) then
2076 c write(*,*) 'Zero K ', t(i,1,1,e,ifld_k-1)
2077  xkey_neg = min(xkey_neg,t(i,1,1,e,ifld_k-1))
2078  nkey_neg = nkey_neg + 1
2079 c write(*,*) 'Neg KEY ', nkey_neg, xkey_neg
2080  t(i,1,1,e,ifld_k-1) = 0.01*abs(t(i,1,1,e,ifld_k-1))
2081  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
2082  endif
2083 
2084  elseif(iflim_omeg.eq.1) then
2085 
2086  if(omega.lt.0.0) then
2087  write(*,*) 'OMEG tot is neg', omega
2088  omega = 0.01*abs(omega)
2089  endif
2090 
2091  if(k.lt.0.0) then
2092  write(*,*) 'K is neg', k
2093  k = 0.01*abs(k)
2094  endif
2095 
2096  endif
2097 
2098  re_t = rho * k /(mu * omega + tiny)
2099  alp_str = alpinf_str * (alp0_str + (re_t/r_k))
2100  $ / (1.+(re_t/r_k))
2101  mu_t = rho * alp_str*k/(omega + tiny) ! should multiply these by rho!!!
2102 
2103  mut(i,1,1,e) = mu_t
2104  mutsk(i,1,1,e) = mu_t / sigma_k
2105  mutso(i,1,1,e) = mu_t / sigma_omega
2106  enddo
2107 
2108  enddo
2109 
2110  if(loglevel.gt.2) then
2111  nome_neg =iglsum(nome_neg,1)
2112  nkey_neg =iglsum(nkey_neg,1)
2113  xome_neg = glmin(xome_neg,1)
2114  xkey_neg = glmin(xkey_neg,1)
2115 
2116  if(nid.eq.0 .and. nome_neg.gt.0)
2117  $ write(*,*) 'Neg Omega ', nome_neg, xome_neg
2118  if(nid.eq.0 .and. nkey_neg.gt.0)
2119  $ write(*,*) 'Neg TKE ', nkey_neg, xkey_neg
2120  endif
2121 
2122  return
2123  end
2124 c-----------------------------------------------------------------------
2125  subroutine rans_komgsst_stndrd_eddy
2126 c
2127 c Compute RANS source terms and diffusivities on an
2128 c element-by-element basis
2129 c
2130  include 'SIZE'
2131  include 'TOTAL'
2132  include 'RANS_KOMG'
2133 
2134  parameter(lxyz=lx1*ly1*lz1)
2135 
2136  real k_x(lxyz),k_y(lxyz),k_z(lxyz)
2137  $ , o_x(lxyz),o_y(lxyz),o_z(lxyz)
2138  $ ,omp_x(lxyz), omp_y(lxyz), omp_z(lxyz)
2139 
2140  common /storesom/ st_mag2(lx1*ly1*lz1,lelv)
2141  $ , om_mag2(lx1*ly1*lz1,lelv)
2142  $ , oiojsk(lx1*ly1*lz1,lelv)
2143  $ , divq(lx1*ly1*lz1,lelv)
2144 
2145  integer e
2146  real k,mu_t,nu_t,mu,nu, kv_min,omeg_max
2147 
2148 c====turbulence constants==========
2149 
2150 c Turbulent viscosity constants
2151  pr_t = coeffs( 1)
2152  sigk1 = coeffs( 2)
2153  sigom1 = coeffs( 3)
2154 
2155 c Low Reynolds number correction constants
2156  alpinf_str = coeffs( 4)
2157  r_k = coeffs( 5)
2158  beta1 = coeffs( 6)
2159  alp0_str = coeffs( 7)
2160 
2161 c Dissipation of K constants
2162  beta_str = coeffs( 8)
2163  gamma1 = coeffs( 9)
2164  r_b = coeffs(10)
2165  akk = coeffs(11)
2166 
2167 c Production of omega constants
2168  alpha_0 = coeffs(12)
2169  r_w = coeffs(13)
2170 
2171 c Dissipation of omega constants
2172 c beta_0 = defined earlier
2173 
2174  kv_min = coeffs(14)
2175  omeg_max = coeffs(15)
2176  tiny = coeffs(16)
2177 
2178 c additional SST and k and epsilon constants
2179  alp1 = coeffs(17)
2180  beta2 = coeffs(18)
2181  sigk2 = coeffs(19)
2182  sigom2 = coeffs(20)
2183  gamma2 = coeffs(21)
2184 c================================
2185 
2186  call comp_stom (st_mag2, om_mag2, oiojsk, divq)
2187 
2188  nome_neg = 0
2189  nkey_neg = 0
2190  xome_neg = 0.
2191  xkey_neg = 0.
2192 
2193  do e=1,nelv
2194 
2195  call gradm11(k_x, k_y, k_z, t(1,1,1,1,ifld_k -1),e)
2196  call gradm11(omp_x,omp_y,omp_z,t(1,1,1,1,ifld_omega-1),e)
2197 
2198  do i=1,lxyz
2199 
2200  rho = param(1) ! vtrans(i,1,1,e,1)
2201  mu = param(2) ! vdiff (i,1,1,e,1)
2202  nu = mu/rho
2203 
2204 c limits for k, omega
2205 
2206 c solve for omega_pert
2207  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
2208  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
2209 
2210  iflim_omeg = 0 ! limit omega^{prime} 1 limit omega_total
2211 
2212  if (iflim_omeg.eq.0) then
2213  if(t(i,1,1,e,ifld_omega-1).lt.0.0) then
2214 c write(*,*) 'Zero OMEG ', t(i,1,1,e,ifld_omega-1)
2215  xome_neg = min(xome_neg,t(i,1,1,e,ifld_omega-1))
2216  nome_neg = nome_neg + 1
2217 c write(*,*) 'Neg OMEG ', nome_neg, xome_neg
2218  t(i,1,1,e,ifld_omega-1) =0.01*abs(t(i,1,1,e,ifld_omega-1))
2219  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
2220  endif
2221 
2222  if(t(i,1,1,e,ifld_k-1).lt.0.0) then
2223 c write(*,*) 'Zero K ', t(i,1,1,e,ifld_k-1)
2224  xkey_neg = min(xkey_neg,t(i,1,1,e,ifld_k-1))
2225  nkey_neg = nkey_neg + 1
2226 c write(*,*) 'Neg KEY ', nkey_neg, xkey_neg
2227  t(i,1,1,e,ifld_k-1) = 0.01*abs(t(i,1,1,e,ifld_k-1))
2228  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
2229  endif
2230 
2231  elseif(iflim_omeg.eq.1) then
2232 
2233  if(omega.lt.0.0) then
2234  write(*,*) 'OMEG tot is neg', omega
2235  omega = 0.01*abs(omega)
2236  endif
2237 
2238  if(k.lt.0.0) then
2239  write(*,*) 'K is neg', k
2240  k = 0.01*abs(k)
2241  endif
2242 
2243  endif
2244 
2245  expn = -2.0
2246  o_x(i)= omp_x(i)+expn * dfdx_omegb(i,1,1,e) *f_omegb(i,1,1,e)
2247  o_y(i)= omp_y(i)+expn * dfdy_omegb(i,1,1,e) *f_omegb(i,1,1,e)
2248  if(if3d)
2249  $ o_z(i)= omp_z(i)+expn * dfdz_omegb(i,1,1,e) *f_omegb(i,1,1,e)
2250 
2251  st_magn = sqrt(st_mag2(i,e))
2252  om_magn = sqrt(om_mag2(i,e))
2253 
2254 ! calculate del k * del omega / omega
2255 
2256  if(if3d)then
2257  xk = (k_x(i)*o_x(i) + k_y(i)*o_y(i) + k_z(i)*o_z(i))
2258  $ / (omega+tiny)
2259  else
2260  xk = (k_x(i)*o_x(i) + k_y(i)*o_y(i))
2261  $ / (omega+tiny)
2262  endif
2263 
2264 ! calculate F2 based on arg2
2265 
2266  yw = ywd(i,1,1,e)
2267  ywm1 = ywdm1(i,1,1,e)
2268  ywm2 = ywm1*ywm1
2269  arg2_1 = sqrt(abs(k)) * ywm1 / omega / beta_str
2270  arg2_2 = 500.0*nu * ywm2 / omega
2271  arg2 = max(2.0*arg2_1, arg2_2)
2272  fun2 = tanh(arg2 * arg2)
2273 
2274 ! calculate F1 based on arg1
2275 
2276  argcd = 2.0 * rho * sigom2 * xk
2277  cdkom = max(argcd, 1.0e-10)
2278  arg1_1 = max( arg2_1, arg2_2)
2279  arg1_2 = 4.0 * rho * sigom2 * k * ywm2 / cdkom
2280  arg1 = min( arg1_1, arg1_2)
2281  fun1 = tanh(arg1 * arg1 * arg1 * arg1)
2282 
2283 ! calculate mu_t
2284  argn_1 = alp1*(omega + tiny)
2285  argn_2 = fun2*st_magn ! this can also be St_magn
2286 
2287  mu_t = 0.0
2288  denom = max(argn_1, argn_2)
2289  if(yw.ne.0) mu_t = rho * alp1 * k / denom
2290 
2291  sigk = fun1 * sigk1 + (1.0 - fun1) * sigk2
2292  sigom = fun1 * sigom1 + (1.0 - fun1) * sigom2
2293 
2294  mut(i,1,1,e) = mu_t
2295  mutsk(i,1,1,e) = mu_t * sigk
2296  mutso(i,1,1,e) = mu_t * sigom
2297 
2298  enddo
2299 
2300  enddo
2301 
2302  if(loglevel.gt.2) then
2303  nome_neg =iglsum(nome_neg,1)
2304  nkey_neg =iglsum(nkey_neg,1)
2305  xome_neg = glmin(xome_neg,1)
2306  xkey_neg = glmin(xkey_neg,1)
2307 
2308  if(nid.eq.0 .and. nome_neg.gt.0)
2309  $ write(*,*) 'Neg Omega ', nome_neg, xome_neg
2310  if(nid.eq.0 .and. nkey_neg.gt.0)
2311  $ write(*,*) 'Neg TKE ', nkey_neg, xkey_neg
2312  endif
2313 
2314  return
2315  end
2316 c-----------------------------------------------------------------------
2317  subroutine rans_komgsst_lowre_eddy
2318 c
2319 c Compute RANS source terms and diffusivities on an
2320 c element-by-element basis
2321 c
2322  include 'SIZE'
2323  include 'TOTAL'
2324  include 'RANS_KOMG'
2325 
2326  parameter(lxyz=lx1*ly1*lz1)
2327 
2328  real k_x(lxyz),k_y(lxyz),k_z(lxyz)
2329  $ , o_x(lxyz),o_y(lxyz),o_z(lxyz)
2330  $ ,omp_x(lxyz), omp_y(lxyz), omp_z(lxyz)
2331 
2332  common /storesom/ st_mag2(lx1*ly1*lz1,lelv)
2333  $ , om_mag2(lx1*ly1*lz1,lelv)
2334  $ , oiojsk(lx1*ly1*lz1,lelv)
2335  $ , divq(lx1*ly1*lz1,lelv)
2336 
2337  integer e
2338  real k,mu_t,nu_t,mu,nu, kv_min,omeg_max
2339 
2340 c====turbulence constants==========
2341 
2342 c Turbulent viscosity constants
2343  pr_t = coeffs( 1)
2344  sigk1 = coeffs( 2)
2345  sigom1 = coeffs( 3)
2346 
2347 c Low Reynolds number correction constants
2348  alpinf_str = coeffs( 4)
2349  r_k = coeffs( 5)
2350  beta1 = coeffs( 6)
2351  alp0_str = coeffs( 7)
2352 
2353 c Dissipation of K constants
2354  beta_str = coeffs( 8)
2355  alp_inf = coeffs( 9)
2356  r_b = coeffs(10)
2357  akk = coeffs(11)
2358 
2359 c Production of omega constants
2360  alpha_0 = coeffs(12)
2361  r_w = coeffs(13)
2362 
2363 c Dissipation of omega constants
2364 c beta_0 = defined earlier
2365 
2366  kv_min = coeffs(14)
2367  omeg_max = coeffs(15)
2368  tiny = coeffs(16)
2369 
2370 c additional SST and k and epsilon constants
2371  alp1 = coeffs(17)
2372  beta2 = coeffs(18)
2373  sigk2 = coeffs(19)
2374  sigom2 = coeffs(20)
2375  gamma2 = coeffs(21)
2376 c================================
2377 
2378  call comp_stom (st_mag2, om_mag2, oiojsk, divq)
2379 
2380  nome_neg = 0
2381  nkey_neg = 0
2382  xome_neg = 0.
2383  xkey_neg = 0.
2384 
2385  do e=1,nelv
2386 
2387  call gradm11(k_x, k_y, k_z, t(1,1,1,1,ifld_k -1),e)
2388  call gradm11(omp_x,omp_y,omp_z,t(1,1,1,1,ifld_omega-1),e)
2389 
2390  do i=1,lxyz
2391 
2392  rho = param(1) ! vtrans(i,1,1,e,1)
2393  mu = param(2) ! vdiff (i,1,1,e,1)
2394  nu = mu/rho
2395 
2396 c limits for k, omega
2397 
2398 c solve for omega_pert
2399  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
2400  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
2401 
2402  iflim_omeg = 0 ! limit omega^{prime} 1 limit omega_total
2403 
2404  if (iflim_omeg.eq.0) then
2405  if(t(i,1,1,e,ifld_omega-1).lt.0.0) then
2406 c write(*,*) 'Zero OMEG ', t(i,1,1,e,ifld_omega-1)
2407  xome_neg = min(xome_neg,t(i,1,1,e,ifld_omega-1))
2408  nome_neg = nome_neg + 1
2409 c write(*,*) 'Neg OMEG ', nome_neg, xome_neg
2410  t(i,1,1,e,ifld_omega-1) =0.01*abs(t(i,1,1,e,ifld_omega-1))
2411  omega = t(i,1,1,e,ifld_omega-1) + f_omegb(i,1,1,e) ! Current k & omega values
2412  endif
2413 
2414  if(t(i,1,1,e,ifld_k-1).lt.0.0) then
2415 c write(*,*) 'Zero K ', t(i,1,1,e,ifld_k-1)
2416  xkey_neg = min(xkey_neg,t(i,1,1,e,ifld_k-1))
2417  nkey_neg = nkey_neg + 1
2418 c write(*,*) 'Neg KEY ', nkey_neg, xkey_neg
2419  t(i,1,1,e,ifld_k-1) = 0.01*abs(t(i,1,1,e,ifld_k-1))
2420  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
2421  endif
2422 
2423  elseif(iflim_omeg.eq.1) then
2424 
2425  if(omega.lt.0.0) then
2426  write(*,*) 'OMEG tot is neg', omega
2427  omega = 0.01*abs(omega)
2428  endif
2429 
2430  if(k.lt.0.0) then
2431  write(*,*) 'K is neg', k
2432  k = 0.01*abs(k)
2433  endif
2434 
2435  endif
2436 
2437  expn = -2.0
2438  o_x(i)= omp_x(i)+expn * dfdx_omegb(i,1,1,e) *f_omegb(i,1,1,e)
2439  o_y(i)= omp_y(i)+expn * dfdy_omegb(i,1,1,e) *f_omegb(i,1,1,e)
2440  if(if3d)
2441  $ o_z(i)= omp_z(i)+expn * dfdz_omegb(i,1,1,e) *f_omegb(i,1,1,e)
2442 
2443  st_magn = sqrt(st_mag2(i,e))
2444  om_magn = sqrt(om_mag2(i,e))
2445 
2446 ! calculate del k * del omega / omega
2447 
2448  if(if3d)then
2449  xk = (k_x(i)*o_x(i) + k_y(i)*o_y(i) + k_z(i)*o_z(i))
2450  $ / (omega+tiny)
2451  else
2452  xk = (k_x(i)*o_x(i) + k_y(i)*o_y(i))
2453  $ / (omega+tiny)
2454  endif
2455 
2456 c ------------ begin low Re part of k-omega -------------------
2457 
2458  re_t = rho * k /(mu * omega + tiny)
2459  alp_str = alpinf_str * (alp0_str + (re_t/r_k))
2460  $ / (1.+(re_t/r_k))
2461 c ------------ end low Re part of k-omega -------------------
2462 
2463 ! calculate F2 based on arg2
2464 
2465  yw = ywd(i,1,1,e)
2466  ywm1 = ywdm1(i,1,1,e)
2467  ywm2 = ywm1*ywm1
2468  arg2_1 = sqrt(abs(k)) * ywm1 / omega / beta_str
2469  arg2_2 = 500.0*nu * ywm2 / omega
2470  arg2 = max(2.0*arg2_1, arg2_2)
2471  fun2 = tanh(arg2 * arg2)
2472 
2473 ! calculate F1 based on arg1
2474 
2475  argcd = 2.0 * rho * sigom2 * xk
2476  cdkom = max(argcd, 1.0e-10)
2477  arg1_1 = max( arg2_1, arg2_2)
2478  arg1_2 = 4.0 * rho * sigom2 * k * ywm2 / cdkom
2479  arg1 = min( arg1_1, arg1_2)
2480  fun1 = tanh(arg1 * arg1 * arg1 * arg1)
2481 
2482 ! calculate mu_t
2483  argn_1 = alp1*(omega + tiny) / alp_str
2484  argn_2 = fun2*st_magn ! this can also be St_magn
2485 
2486  mu_t = 0.0
2487  denom = max(argn_1, argn_2)
2488  if(yw.ne.0) mu_t = rho * alp1 * k / denom
2489 
2490  sigk = fun1 * sigk1 + (1.0 - fun1) * sigk2
2491  sigom = fun1 * sigom1 + (1.0 - fun1) * sigom2
2492 
2493  mut(i,1,1,e) = mu_t
2494  mutsk(i,1,1,e) = mu_t * sigk
2495  mutso(i,1,1,e) = mu_t * sigom
2496 
2497  enddo
2498 
2499  enddo
2500 
2501  if(loglevel.gt.2) then
2502  nome_neg =iglsum(nome_neg,1)
2503  nkey_neg =iglsum(nkey_neg,1)
2504  xome_neg = glmin(xome_neg,1)
2505  xkey_neg = glmin(xkey_neg,1)
2506 
2507  if(nid.eq.0 .and. nome_neg.gt.0)
2508  $ write(*,*) 'Neg Omega ', nome_neg, xome_neg
2509  if(nid.eq.0 .and. nkey_neg.gt.0)
2510  $ write(*,*) 'Neg TKE ', nkey_neg, xkey_neg
2511  endif
2512 
2513  return
2514  end
2515 c-----------------------------------------------------------------------
2516  subroutine rans_komg_stndrd_eddy_noreg
2517 c
2518 c Compute RANS source terms and diffusivities on an
2519 c element-by-element basis
2520 c
2521  include 'SIZE'
2522  include 'TOTAL'
2523  include 'RANS_KOMG'
2524 
2525  parameter(lxyz=lx1*ly1*lz1)
2526 
2527  integer e
2528  real k,mu_t,nu_t,mu,nu, kv_min,omeg_max
2529 
2530 c Turbulent viscosity constants
2531  pr_t = coeffs( 1)
2532  sigma_k = coeffs( 2)
2533  sigma_omega = coeffs( 3)
2534 
2535 c Low Reynolds number correction constants
2536  alpinf_str = coeffs( 4)
2537  r_k = coeffs( 5)
2538  beta_0 = coeffs( 6)
2539  alp0_str = coeffs( 7)
2540 
2541 c Dissipation of K constants
2542  betainf_str = coeffs( 8)
2543  alp_inf = coeffs( 9)
2544  r_b = coeffs(10)
2545  akk = coeffs(11)
2546 
2547 c Production of omega constants
2548  alpha_0 = coeffs(12)
2549  r_w = coeffs(13)
2550 
2551 c Dissipation of omega constants
2552 c beta_0 = defined earlier
2553 
2554  kv_min = coeffs(14)
2555  omeg_max = coeffs(15)
2556  tiny = coeffs(16)
2557 c================================
2558 
2559  nome_neg = 0
2560  nkey_neg = 0
2561  xome_neg = 0.
2562  xkey_neg = 0.
2563 
2564  do e=1,nelv
2565  do i=1,lxyz
2566 
2567  rho = param(1) ! vtrans(i,1,1,e,1)
2568  mu = param(2) ! vdiff (i,1,1,e,1)
2569  nu = mu/rho
2570 
2571 c limits for k, omega
2572 
2573 c solve for omega_pert
2574  omega = t(i,1,1,e,ifld_omega-1) ! Current k & omega values
2575  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
2576 
2577  iflim_omeg = 0 ! limit omega^{prime} 1 limit omega_total
2578 
2579  if (iflim_omeg.eq.0) then
2580  if(t(i,1,1,e,ifld_omega-1).lt.0.0) then
2581 c write(*,*) 'Zero OMEG ', t(i,1,1,e,ifld_omega-1)
2582  xome_neg = min(xome_neg,t(i,1,1,e,ifld_omega-1))
2583  nome_neg = nome_neg + 1
2584 c write(*,*) 'Neg OMEG ', nome_neg, xome_neg
2585  t(i,1,1,e,ifld_omega-1) =0.01*abs(t(i,1,1,e,ifld_omega-1))
2586  omega = t(i,1,1,e,ifld_omega-1) ! Current k & omega values
2587  endif
2588 
2589  if(t(i,1,1,e,ifld_k-1).lt.0.0) then
2590 c write(*,*) 'Zero K ', t(i,1,1,e,ifld_k-1)
2591  xkey_neg = min(xkey_neg,t(i,1,1,e,ifld_k-1))
2592  nkey_neg = nkey_neg + 1
2593 c write(*,*) 'Neg KEY ', nkey_neg, xkey_neg
2594  t(i,1,1,e,ifld_k-1) = 0.01*abs(t(i,1,1,e,ifld_k-1))
2595  k = t(i,1,1,e,ifld_k -1) ! from previous timestep
2596  endif
2597 
2598  elseif(iflim_omeg.eq.1) then
2599 
2600  if(omega.lt.0.0) then
2601  write(*,*) 'OMEG tot is neg', omega
2602  omega = 0.01*abs(omega)
2603  endif
2604 
2605  if(k.lt.0.0) then
2606  write(*,*) 'K is neg', k
2607  k = 0.01*abs(k)
2608  endif
2609 
2610  endif
2611 
2612  alp_str = alpinf_str
2613  mu_t = rho * alp_str*k/(omega + tiny) ! should multiply these by rho!!!
2614 
2615  mut(i,1,1,e) = mu_t
2616  mutsk(i,1,1,e) = mu_t / sigma_k
2617  mutso(i,1,1,e) = mu_t / sigma_omega
2618  enddo
2619 
2620  enddo
2621 
2622  if(loglevel.gt.2) then
2623  nome_neg =iglsum(nome_neg,1)
2624  nkey_neg =iglsum(nkey_neg,1)
2625  xome_neg = glmin(xome_neg,1)
2626  xkey_neg = glmin(xkey_neg,1)
2627 
2628  if(nid.eq.0 .and. nome_neg.gt.0)
2629  $ write(*,*) 'Neg Omega ', nome_neg, xome_neg
2630  if(nid.eq.0 .and. nkey_neg.gt.0)
2631  $ write(*,*) 'Neg TKE ', nkey_neg, xkey_neg
2632  endif
2633 
2634  return
2635  end
2636 c-----------------------------------------------------------------------
2637  subroutine comp_stom(St_mag2,Om_mag2,OiOjSk,DivQ)
2638 c
2639 c Compute the square of the magnitude of the stress and rotation tensors
2640 c St_mag2=2Sij*Sij=S'ij*S'ij/2, S'ij=2Sij, Sij=(dui/dxj+duj/dxi)/2
2641 c Om_mag2=2Oij*Oij=O'ij*O'ij/2, O'ij=OSij, Oij=(dui/dxj-duj/dxi)/2
2642 c
2643  include 'SIZE'
2644  include 'TOTAL'
2645 c
2646  integer e
2647 
2648  parameter(lt =lx1*ly1*lz1*lelv)
2649  parameter(lxyz=lx1*ly1*lz1)
2650 
2651  common /scruz/ sij(lx1*ly1*lz1,6,lelv)
2652  $ , oij(lx1*ly1*lz1,lelv,3)
2653 
2654  real St_mag2(lx1*ly1*lz1,lelv)
2655  $ , om_mag2(lx1*ly1*lz1,lelv)
2656  $ , oiojsk(lx1*ly1*lz1,lelv)
2657  $ , divq(lx1*ly1*lz1,lelv)
2658 
2659  logical iflmc, ifdss
2660 c
2661  thqrt = 0.75
2662 
2663  iflmc = .false.! to add lowmach correction to 2Sij (2Sij - 2Q/3*delta_ij)
2664  ifdss = .true. ! to dssum sij and oij
2665  nij = 3
2666  if (if3d.or.ifaxis) nij=6
2667 
2668  if(nid.eq.0 .and. loglevel.gt.2) write(*,*) 'updating StOm '
2669  call comp_sij_oij (sij, oij, nij, vx, vy, vz, iflmc, ifdss) ! S'ij=2Sij, O'ij=2Oij
2670  call comp_sij_oij_mag2(sij, oij, nij, st_mag2, om_mag2, oiojsk) ! St_mag2=2S^2=S'^2/2
2671 c ! Om_mag2=2O^2=O'^2/2
2672  do e = 1, nelv
2673  call add3 (divq(1,e), sij(1,1,e), sij(1,2,e), lxyz)
2674  if(if3d.or.ifaxis)
2675  $ call add2 (divq(1,e), sij(1,3,e), lxyz)
2676  if(iflmc) call cmult(divq(1,e), thqrt, lxyz)
2677  enddo
2678 
2679  return
2680  end
2681 c-----------------------------------------------------------------------
2682  subroutine comp_sij_oij(sij,oij,nij,u,v,w,iflmc,ifdss)
2683 c du_i du_j
2684 c Compute the stress tensor S_ij := ---- + ----
2685 c du_j du_i
2686 c
2687  include 'SIZE'
2688  include 'TOTAL'
2689 c
2690  integer e
2691  logical iflmc, ifdss
2692 c
2693  parameter(lt =lx1*ly1*lz1*lelv)
2694  parameter(lxyz=lx1*ly1*lz1)
2695 
2696  real sij(lx1*ly1*lz1,nij,lelv), oij(lx1*ly1*lz1,lelv,3)
2697  real u (lx1*ly1*lz1,lelv)
2698  real v (lx1*ly1*lz1,lelv)
2699  real w (lx1*ly1*lz1,lelv)
2700 
2701  real ur(lxyz),us(lxyz),ut(lxyz)
2702  $ ,vr(lxyz),vs(lxyz),vt(lxyz)
2703  $ ,wr(lxyz),ws(lxyz),wt(lxyz)
2704 
2705  real j ! Inverse Jacobian
2706 
2707  n = lx1-1 ! Polynomial degree
2708  nxyz = lx1*ly1*lz1
2709 
2710  if (if3d) then ! 3D CASE
2711  do e=1,nelv
2712  call local_grad3(ur,us,ut,u,n,e,dxm1,dxtm1)
2713  call local_grad3(vr,vs,vt,v,n,e,dxm1,dxtm1)
2714  call local_grad3(wr,ws,wt,w,n,e,dxm1,dxtm1)
2715 
2716  do i=1,nxyz
2717 
2718  j = jacmi(i,e)
2719 
2720  sij(i,1,e) = j* ! du/dx + du/dx
2721  $ 2*(ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e)+ut(i)*txm1(i,1,1,e))
2722 
2723  sij(i,2,e) = j* ! dv/dy + dv/dy
2724  $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e)+vt(i)*tym1(i,1,1,e))
2725 
2726  sij(i,3,e) = j* ! dw/dz + dw/dz
2727  $ 2*(wr(i)*rzm1(i,1,1,e)+ws(i)*szm1(i,1,1,e)+wt(i)*tzm1(i,1,1,e))
2728 
2729  sij(i,4,e) = j* ! du/dy + dv/dx
2730  $ (ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e)+ut(i)*tym1(i,1,1,e) +
2731  $ vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e)+vt(i)*txm1(i,1,1,e) )
2732  oij(i,e,3) = j* ! dv/dx - du/dy
2733  $ (vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e)+vt(i)*txm1(i,1,1,e) -
2734  $ ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e)+ut(i)*tym1(i,1,1,e) )
2735 
2736  sij(i,5,e) = j* ! dv/dz + dw/dy
2737  $ (wr(i)*rym1(i,1,1,e)+ws(i)*sym1(i,1,1,e)+wt(i)*tym1(i,1,1,e) +
2738  $ vr(i)*rzm1(i,1,1,e)+vs(i)*szm1(i,1,1,e)+vt(i)*tzm1(i,1,1,e) )
2739  oij(i,e,1) = j* ! dw/dy - dv/dz
2740  $ (wr(i)*rym1(i,1,1,e)+ws(i)*sym1(i,1,1,e)+wt(i)*tym1(i,1,1,e) -
2741  $ vr(i)*rzm1(i,1,1,e)+vs(i)*szm1(i,1,1,e)+vt(i)*tzm1(i,1,1,e) )
2742 
2743  sij(i,6,e) = j* ! du/dz + dw/dx
2744  $ (ur(i)*rzm1(i,1,1,e)+us(i)*szm1(i,1,1,e)+ut(i)*tzm1(i,1,1,e) +
2745  $ wr(i)*rxm1(i,1,1,e)+ws(i)*sxm1(i,1,1,e)+wt(i)*txm1(i,1,1,e) )
2746  oij(i,e,2) = j* ! du/dz - dw/dx
2747  $ (ur(i)*rzm1(i,1,1,e)+us(i)*szm1(i,1,1,e)+ut(i)*tzm1(i,1,1,e) -
2748  $ wr(i)*rxm1(i,1,1,e)+ws(i)*sxm1(i,1,1,e)+wt(i)*txm1(i,1,1,e) )
2749 
2750 
2751  enddo
2752  enddo
2753 
2754  elseif (ifaxis) then ! AXISYMMETRIC CASE
2755 
2756 c
2757 c Notation: ( 2 x Acheson, p. 353)
2758 c Cylindrical
2759 c Nek5k Coordinates
2760 c
2761 c x z
2762 c y r
2763 c z theta
2764 c
2765 
2766  do e=1,nelv
2767  call setaxdy ( ifrzer(e) ) ! change dytm1 if on-axis
2768  call local_grad2(ur,us,u,n,e,dxm1,dytm1)
2769  call local_grad2(vr,vs,v,n,e,dxm1,dytm1)
2770  call local_grad2(wr,ws,w,n,e,dxm1,dytm1)
2771 
2772  do i=1,nxyz
2773  j = jacmi(i,e)
2774  r = ym1(i,1,1,e) ! Cyl. Coord:
2775 
2776  sij(i,1,e) = j* ! du/dx + du/dx ! e_zz
2777  $ 2*(ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e))
2778 
2779  sij(i,2,e) = j* ! dv/dy + dv/dy ! e_rr
2780  $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
2781 
2782  if (r.gt.0) then ! e_@@
2783  sij(i,3,e) = 2.*v(i,e)/r ! v / r ! corrected AT: factor of 2, 10/30/18
2784  else
2785  sij(i,3,e) = j* ! L'Hopital's rule: e_@@ = 2dv/dr
2786  $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
2787  endif
2788 
2789  sij(i,4,e) = j* ! du/dy + dv/dx ! e_zr
2790  $ ( ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e) +
2791  $ vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e) )
2792  oij(i,e,3) = j* ! dv/dx - du/dy
2793  $ ( vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e) -
2794  $ ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e) )
2795 
2796  if (r.gt.0) then ! e_r@
2797  sij(i,5,e) = j* ! dw/dy
2798  $ ( wr(i)*rym1(i,1,1,e)+ws(i)*sym1(i,1,1,e) )
2799  $ - w(i,e) / r
2800  oij(i,e,1) = j* ! dw/dy
2801  $ ( wr(i)*rym1(i,1,1,e)+ws(i)*sym1(i,1,1,e) )
2802  $ + w(i,e) / r
2803  else
2804  sij(i,5,e) = 0
2805  oij(i,e,2) = j* ! L'Hopital's rule: e_r@ = 2dw/dr
2806  $ 2*(wr(i)*rym1(i,1,1,e)+ws(i)*sym1(i,1,1,e))
2807  endif
2808 
2809  sij(i,6,e) = j* ! dw/dx ! e_@z
2810  $ ( wr(i)*rxm1(i,1,1,e)+ws(i)*sxm1(i,1,1,e) )
2811  oij(i,e,2) = j* !-dw/dx
2812  $ (-wr(i)*rxm1(i,1,1,e)-ws(i)*sxm1(i,1,1,e) )
2813 
2814  enddo
2815  enddo
2816 
2817  else ! 2D CASE
2818 
2819  do e=1,nelv
2820  call local_grad2(ur,us,u,n,e,dxm1,dxtm1)
2821  call local_grad2(vr,vs,v,n,e,dxm1,dxtm1)
2822 
2823  do i=1,nxyz
2824  j = jacmi(i,e)
2825 
2826  sij(i,1,e) = j* ! du/dx + du/dx
2827  $ 2*(ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e))
2828  oij(i,e,2) = 0.
2829 
2830  sij(i,2,e) = j* ! dv/dy + dv/dy
2831  $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
2832  oij(i,e,3) = 0.
2833 
2834  sij(i,3,e) = j* ! du/dy + dv/dx
2835  $ (ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e) +
2836  $ vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e) )
2837  oij(i,e,1) = j* ! dv/dx - du/dy
2838  $ (vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e) -
2839  $ ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e) )
2840 
2841  enddo
2842  enddo
2843  endif
2844 
2845  if(iflomach .and. iflmc) then
2846 
2847  onethird = 1./3.
2848  if(if3d .or. ifaxis) then
2849  do e=1,nelv
2850  do i=1,nxyz
2851  trs = sij(i,1,e) + sij(i,2,e) + sij(i,3,e) ! 2(du/dx + dv/dy + dw/dz or v/r) in axisym
2852  sij(i,1,e) = sij(i,1,e) - onethird*trs ! 2S - (2/3)Q = S'-(1/3)tr(S')
2853  enddo
2854  enddo
2855  else
2856  do e=1,nelv
2857  do i=1,nxyz
2858  trs = sij(i,1,e) + sij(i,2,e) ! 2(du/dx + dv/dy)
2859  sij(i,1,e) = sij(i,1,e) - onethird*trs ! 2S - (2/3)Q = S'-(1/3)tr(S')
2860  enddo
2861  enddo
2862  endif
2863 
2864  endif
2865 
2866  ifielt = ifield
2867  ifield = 1
2868  if(ifdss) call dssum_sij_oij(sij, oij, nij)
2869  ifield = ifielt
2870 
2871  return
2872  end
2873 c-----------------------------------------------------------------------
2874  subroutine dssum_sij_oij(sij,oij,nij)
2875 c
2876  include 'SIZE'
2877  include 'TOTAL'
2878 c
2879  integer e
2880 c
2881  parameter(lt =lx1*ly1*lz1*lelv)
2882  parameter(lxyz=lx1*ly1*lz1)
2883 
2884  real sij(lx1*ly1*lz1,nij,lelv), oij(lx1*ly1*lz1,lelv,3)
2885 
2886  common /scrsij/ work1(lx1*ly1*lz1,lelv), work2(lx1*ly1*lz1,lelv)
2887  $ , work3(lx1*ly1*lz1,lelv)
2888 
2889  nxyz = lx1*ly1*lz1
2890  ntot = nxyz*nelv
2891 
2892  if (if3d .or. ifaxis) then
2893  call opcolv (oij(1,1,1), oij(1,1,2), oij(1,1,3), bm1)
2894  call opdssum (oij(1,1,1), oij(1,1,2), oij(1,1,3) )
2895  call opcolv (oij(1,1,1), oij(1,1,2), oij(1,1,3),binvm1)
2896  do e=1,nelv
2897  call copy (work1(1,e), sij(1,1,e),nxyz)
2898  call copy (work2(1,e), sij(1,2,e),nxyz)
2899  call copy (work3(1,e), sij(1,3,e),nxyz)
2900  enddo
2901  call opcolv (work1, work2, work3, bm1)
2902  call opdssum (work1, work2, work3 )
2903  call opcolv (work1, work2, work3,binvm1)
2904  do e=1,nelv
2905  call copy (sij(1,1,e),work1(1,e), nxyz)
2906  call copy (sij(1,2,e),work2(1,e), nxyz)
2907  call copy (sij(1,3,e),work3(1,e), nxyz)
2908  enddo
2909  do e=1,nelv
2910  call copy (work1(1,e), sij(1,4,e),nxyz)
2911  call copy (work2(1,e), sij(1,5,e),nxyz)
2912  call copy (work3(1,e), sij(1,6,e),nxyz)
2913  enddo
2914  call opcolv (work1, work2, work3, bm1)
2915  call opdssum (work1, work2, work3 )
2916  call opcolv (work1, work2, work3,binvm1)
2917  do e=1,nelv
2918  call copy (sij(1,4,e),work1(1,e), nxyz)
2919  call copy (sij(1,5,e),work2(1,e), nxyz)
2920  call copy (sij(1,6,e),work3(1,e), nxyz)
2921  enddo
2922  else
2923  call col2 (oij(1,1,1),bm1 ,ntot)
2924  call dssum (oij(1,1,1),lx1,ly1,lz1)
2925  call col2 (oij(1,1,1),binvm1,ntot)
2926  do e=1,nelv
2927  call copy (work1(1,e), sij(1,1,e),nxyz)
2928  call copy (work2(1,e), sij(1,2,e),nxyz)
2929  call copy (work3(1,e), sij(1,3,e),nxyz)
2930  enddo
2931  call opcolv (work1, work2, work3, bm1)
2932  call opdssum (work1, work2, work3 )
2933  call opcolv (work1, work2, work3,binvm1)
2934  do e=1,nelv
2935  call copy (sij(1,1,e),work1(1,e), nxyz)
2936  call copy (sij(1,2,e),work2(1,e), nxyz)
2937  call copy (sij(1,3,e),work3(1,e), nxyz)
2938  enddo
2939  endif
2940 
2941  return
2942  end
2943 c-----------------------------------------------------------------------
2944  subroutine comp_sij_oij_mag2(sij,oij,nij,St_mag2,Om_mag2,OiOjSk)
2945 c
2946 c Compute the square of the magnitude of the stress and rotation tensors
2947 c St_mag2=2Sij*Sij=S'ij*S'ij/2, S'ij=2Sij, Sij=(dui/dxj+duj/dxi)/2
2948 c Om_mag2=2Oij*Oij=O'ij*O'ij/2, O'ij=OSij, Oij=(dui/dxj-duj/dxi)/2
2949 c
2950  include 'SIZE'
2951  include 'TOTAL'
2952 c
2953  integer e
2954 c
2955  parameter(lt =lx1*ly1*lz1*lelv)
2956  parameter(lxyz=lx1*ly1*lz1)
2957 
2958  real sij(lx1*ly1*lz1,nij,lelv), oij(lx1*ly1*lz1,lelv,3)
2959 
2960  real St_mag2(lx1*ly1*lz1,lelv)
2961  $ , om_mag2(lx1*ly1*lz1,lelv)
2962  $ , oiojsk(lx1*ly1*lz1,lelv)
2963 
2964  common /scrsij/ work1(lx1*ly1*lz1,lelv), work2(lx1*ly1*lz1,lelv)
2965  $ , work3(lx1*ly1*lz1,lelv)
2966  real tmp1(lxyz), tmp2(lxyz), tmp3(lxyz)
2967 
2968  nxyz = lx1*ly1*lz1
2969  ntot = nxyz*nelv
2970  two = 2.
2971  onehalf = 0.5
2972  onethird= 1./3.
2973  oneeight= 1./8.
2974  beta = onehalf-onethird
2975 
2976  if (if3d .or. ifaxis) then ! 3D CASE or axisymmetric
2977 
2978  do e=1,nelv
2979  call col3 (st_mag2(1,e), sij(1,1,e), sij(1,1,e),nxyz)
2980  call addcol3 (st_mag2(1,e), sij(1,2,e), sij(1,2,e),nxyz)
2981  call addcol3 (st_mag2(1,e), sij(1,3,e), sij(1,3,e),nxyz)
2982  call col3 (work1(1,e), sij(1,4,e), sij(1,4,e),nxyz)
2983  call addcol3 (work1(1,e), sij(1,5,e), sij(1,5,e),nxyz)
2984  call addcol3 (work1(1,e), sij(1,6,e), sij(1,6,e),nxyz)
2985  call add2s2 (st_mag2(1,e), work1(1,e), two, nxyz)
2986 c
2987  call col3 (work1(1,e), oij(1,e,1), oij(1,e,1),nxyz)
2988  call col3 (work2(1,e), oij(1,e,2), oij(1,e,2),nxyz)
2989  call col3 (work3(1,e), oij(1,e,3), oij(1,e,3),nxyz)
2990  call add4 (om_mag2(1,e), work1(1,e), work2(1,e)
2991  $ , work3(1,e),nxyz)
2992 
2993  call add3 (tmp1, work2(1,e), work3(1,e), nxyz)
2994  call add3 (tmp2, work1(1,e), work3(1,e), nxyz)
2995  call add3 (tmp3, work1(1,e), work2(1,e), nxyz)
2996  call col3 (oiojsk(1,e),tmp1, sij(1,1,e), nxyz)
2997  call addcol3 (oiojsk(1,e),tmp2, sij(1,2,e), nxyz)
2998  call addcol3 (oiojsk(1,e),tmp3, sij(1,3,e), nxyz)
2999 
3000  call col4 (tmp1, oij(1,e,1), oij(1,e,2), sij(1,4,e), nxyz)
3001  call addcol4 (tmp1, oij(1,e,2), oij(1,e,3), sij(1,5,e), nxyz)
3002  call subcol4 (tmp1, oij(1,e,1), oij(1,e,3), sij(1,6,e), nxyz)
3003  call add2s2 (oiojsk(1,e), tmp1, two, nxyz)
3004 
3005  enddo
3006 
3007  else ! 2D CASE
3008 
3009  do e=1,nelv
3010  call col3 (st_mag2(1,e), sij(1,1,e), sij(1,1,e),nxyz)
3011  call addcol3 (st_mag2(1,e), sij(1,2,e), sij(1,2,e),nxyz)
3012  call col3 (work1(1,e), sij(1,3,e), sij(1,3,e),nxyz)
3013  call add2s2 (st_mag2(1,e), work1(1,e), two, nxyz)
3014 c
3015  call col3 (om_mag2(1,e), oij(1,e,1), oij(1,e,1),nxyz)
3016  call rzero (oiojsk(1,e), nxyz)
3017  enddo
3018 
3019  endif
3020 
3021  call cmult (st_mag2, onehalf, ntot) ! St_mag2=2*Sij*Sij=S'ij*S'ij/2
3022 
3023  return
3024  end
3025 c-----------------------------------------------------------------------
exitt
void exitt()
Definition: comm_mpi.f:604
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.f:535
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
glmin
function glmin(a, n)
Definition: math.f:973
invcol2
subroutine invcol2(a, b, n)
Definition: math.f:73
addcol3
subroutine addcol3(a, b, c, n)
Definition: math.f:654
admcol3
subroutine admcol3(a, b, c, d, n)
Definition: math.f:1556
iglsum
function iglsum(a, n)
Definition: math.f:926
add2
subroutine add2(a, b, n)
Definition: math.f:622
col2
subroutine col2(a, b, n)
Definition: math.f:564
add2s2
subroutine add2s2(a, b, c1, n)
Definition: math.f:690
copy
subroutine copy(a, b, n)
Definition: math.f:260
addcol4
subroutine addcol4(a, b, c, d, n)
Definition: math.f:142
add3
subroutine add3(a, b, c, n)
Definition: math.f:644
col4
subroutine col4(a, b, c, d, n)
Definition: math.f:120
add4
subroutine add4(a, b, c, d, n)
Definition: math.f:728
subcol3
subroutine subcol3(a, b, c, n)
Definition: math.f:186
subcol4
subroutine subcol4(a, b, c, d, n)
Definition: math.f:197
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
rone
subroutine rone(a, n)
Definition: math.f:230
rzero
subroutine rzero(a, n)
Definition: math.f:208
opdssum
subroutine opdssum(a, b, c)
Definition: navier1.f:2582
opcolv
subroutine opcolv(a1, a2, a3, c)
Definition: navier1.f:2418
gradm11
subroutine gradm11(ux, uy, uz, u, e)
Definition: navier2.f:176
cheap_dist
subroutine cheap_dist(d, ifld, b)
Definition: navier5.f:2970
local_grad2
subroutine local_grad2(ur, us, u, N, e, D, Dt)
Definition: navier5.f:448
local_grad3
subroutine local_grad3(ur, us, ut, u, N, e, D, Dt)
Definition: navier5.f:429
distf
subroutine distf(d, ifld, b, dmin, emin, xn, yn, zn)
Definition: navier5.f:3067
gradm1
subroutine gradm1(ux, uy, uz, u)
Definition: navier5.f:463
rans_komgsst_lowre_compute
subroutine rans_komgsst_lowre_compute
Definition: rans_komg.f:981
rans_komg_ksrc
real function rans_komg_ksrc(ix, iy, iz, iel)
Definition: rans_komg.f:59
rans_komg_lowre_eddy
subroutine rans_komg_lowre_eddy
Definition: rans_komg.f:2002
rans_komg_lowre_compute
subroutine rans_komg_lowre_compute
Definition: rans_komg.f:390
rans_komg_init
subroutine rans_komg_init(ifld_k_in, ifld_omega_in, ifcoeffs, coeffs_in, wall_id, ywd_in, model_id)
Definition: rans_komg.f:1558
rans_komg_mutsk
real function rans_komg_mutsk(ix, iy, iz, iel)
Definition: rans_komg.f:25
rans_komg_stndrd_compute
subroutine rans_komg_stndrd_compute
Definition: rans_komg.f:109
rans_komg_mut
real function rans_komg_mut(ix, iy, iz, iel)
Definition: rans_komg.f:2
rans_komgsst_lowre_eddy
subroutine rans_komgsst_lowre_eddy
Definition: rans_komg.f:2318
rans_komgsst_stndrd_compute
subroutine rans_komgsst_stndrd_compute
Definition: rans_komg.f:677
rans_komg_getcoeffs
subroutine rans_komg_getcoeffs(coeffs_in)
Definition: rans_komg.f:45
rans_komg_omgsrc
real function rans_komg_omgsrc(ix, iy, iz, iel)
Definition: rans_komg.f:84
comp_sij_oij_mag2
subroutine comp_sij_oij_mag2(sij, oij, nij, St_mag2, Om_mag2, OiOjSk)
Definition: rans_komg.f:2945
rans_komg_stndrd_eddy
subroutine rans_komg_stndrd_eddy
Definition: rans_komg.f:1880
comp_sij_oij
subroutine comp_sij_oij(sij, oij, nij, u, v, w, iflmc, ifdss)
Definition: rans_komg.f:2683
rans_komg_mutso
real function rans_komg_mutso(ix, iy, iz, iel)
Definition: rans_komg.f:35
rans_komgsst_stndrd_eddy
subroutine rans_komgsst_stndrd_eddy
Definition: rans_komg.f:2126
dssum_sij_oij
subroutine dssum_sij_oij(sij, oij, nij)
Definition: rans_komg.f:2875
comp_stom
subroutine comp_stom(St_mag2, Om_mag2, OiOjSk, DivQ)
Definition: rans_komg.f:2638
rans_komgsst_set_defaultcoeffs
subroutine rans_komgsst_set_defaultcoeffs
Definition: rans_komg.f:1803
rans_komg_set_defaultcoeffs
subroutine rans_komg_set_defaultcoeffs
Definition: rans_komg.f:1740
rans_komg_omegabase
subroutine rans_komg_omegabase
Definition: rans_komg.f:1657
rans_komg_stndrd_compute_noreg
subroutine rans_komg_stndrd_compute_noreg
Definition: rans_komg.f:1330
rans_komg_stndrd_eddy_noreg
subroutine rans_komg_stndrd_eddy_noreg
Definition: rans_komg.f:2517
setaxdy
subroutine setaxdy(ifaxdy)
Definition: subs1.f:2342