Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLO_RoeDissipSecond.F90
Go to the documentation of this file.
1 ! *********************************************************************
2 ! * Rocstar Simulation Suite *
3 ! * Copyright@2015, Illinois Rocstar LLC. All rights reserved. *
4 ! * *
5 ! * Illinois Rocstar LLC *
6 ! * Champaign, IL *
7 ! * www.illinoisrocstar.com *
8 ! * sales@illinoisrocstar.com *
9 ! * *
10 ! * License: See LICENSE file in top level of distribution package or *
11 ! * http://opensource.org/licenses/NCSA *
12 ! *********************************************************************
13 ! *********************************************************************
14 ! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
15 ! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
16 ! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
17 ! * NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
18 ! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
19 ! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
20 ! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
21 ! * USE OR OTHER DEALINGS WITH THE SOFTWARE. *
22 ! *********************************************************************
23 !******************************************************************************
24 !
25 ! Purpose: compute numerical dissipation based on 2nd-order
26 ! Roe`s upwind scheme.
27 !
28 ! Description: none.
29 !
30 ! Input: region = data of current region.
31 !
32 ! Output: region%levels%mixt%diss = dissipative fluxes.
33 !
34 ! Notes: uses MUSCL scheme with kappa=1/3.
35 !
36 !******************************************************************************
37 !
38 ! $Id: RFLO_RoeDissipSecond.F90,v 1.3 2008/12/06 08:44:27 mtcampbe Exp $
39 !
40 ! Copyright: (c) 2003 by the University of Illinois
41 !
42 !******************************************************************************
43 
44 SUBROUTINE rflo_roedissipsecond( region )
45 
46  USE moddatatypes
47  USE moddatastruct, ONLY : t_region
48  USE modinterfaces, ONLY : rflo_getdimensphys, rflo_getcelloffset, &
49  rflo_getnodeoffset
50  USE moderror
51  USE modparameters
52  IMPLICIT NONE
53 
54 #include "Indexing.h"
55 
56 ! ... parameters
57  TYPE(t_region) :: region
58 
59 ! ... loop variables
60  INTEGER :: i, j, k
61 
62 ! ... local variables
63  INTEGER :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
64  INTEGER :: icoff, ijcoff, inoff, ijnoff, ijkc0, ijkcm1, ijkcm2, ijkcp1, ijkn
65  INTEGER :: ilev, indcp, indmol, indsvel
66 
67  REAL(RFREAL) :: beta5, ds, nx, ny, nz, svel, rgas, gaml, gamr, ggm1, gam1, &
68  rl, ul, vl, wl, pl, hl, rr, ur, vr, wr, pr, hr, rav, dd, &
69  dd1, uav, vav, wav, hav, q2a, c2a, cav, uvw, du, eabs1, &
70  eabs2, eabs5, h1, h2, h3, h4, h5, epsentr, delta, limfac, &
71  limfac3, rvolref, eps2(3), eps2n, vola, dvar(5), dvarm(5), &
72  dvarp(5), deltl(5), deltr(5), fd(5)
73  REAL(RFREAL), POINTER :: cv(:,:), dv(:,:), gv(:,:), diss(:,:), &
74  vol(:), si(:,:), sj(:,:), sk(:,:), &
75  sivel(:), sjvel(:), skvel(:)
76 
77 ! ... functions
78  REAL(RFREAL) :: muscl3, af, bf, eps
79 
80 !******************************************************************************
81 ! limiter function
82 
83  muscl3(af,bf,eps) = (bf*(2._rfreal*af*af+eps)+af*(bf*bf+2._rfreal*eps))/ &
84  (2._rfreal*af*af+2._rfreal*bf*bf-af*bf+ &
85  3._rfreal*eps+1.e-30_rfreal)
86 
87  CALL registerfunction( region%global,'RFLO_RoeDissipSecond',&
88  'RFLO_RoeDissipSecond.F90' )
89 
90 ! get dimensions and pointers -------------------------------------------------
91 
92  ilev = region%currLevel
93 
94  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
95  jpcbeg,jpcend,kpcbeg,kpcend )
96  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
97  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
98 
99  cv => region%levels(ilev)%mixt%cv
100  dv => region%levels(ilev)%mixt%dv
101  gv => region%levels(ilev)%mixt%gv
102  vol => region%levels(ilev)%grid%vol
103  si => region%levels(ilev)%grid%si
104  sj => region%levels(ilev)%grid%sj
105  sk => region%levels(ilev)%grid%sk
106  sivel => region%levels(ilev)%grid%siVel
107  sjvel => region%levels(ilev)%grid%sjVel
108  skvel => region%levels(ilev)%grid%skVel
109  diss => region%levels(ilev)%mixt%diss
110 
111 ! get coefficients ------------------------------------------------------------
112 
113  beta5 = 0.5_rfreal*region%mixtInput%betrk(region%irkStep)
114  epsentr = region%mixtInput%epsEntr
115  limfac = region%mixtInput%limFac
116  indsvel = region%levels(ilev)%grid%indSvel
117  indcp = region%levels(ilev)%mixt%indCp
118  indmol = region%levels(ilev)%mixt%indMol
119 
120 ! normalise epsilon^2 for all limited variables (rho, u, v, w, p) -------------
121 
122  limfac3 = limfac*limfac*limfac
123  rvolref = 1._rfreal/region%global%limVolRef**1.5_rfreal
124  eps2(1) = limfac3*region%global%limRef(1)*region%global%limRef(1)*rvolref
125  eps2(2) = limfac3*region%global%limRef(2)*region%global%limRef(2)*rvolref
126  eps2(3) = limfac3*region%global%limRef(3)*region%global%limRef(3)*rvolref
127 
128 ! dissipation in i-direction --------------------------------------------------
129 
130  DO k=kpcbeg,kpcend
131  DO j=jpcbeg,jpcend
132  DO i=ipcbeg,ipcend+1
133  ijkc0 = indijk(i ,j,k,icoff,ijcoff)
134  ijkcm1 = indijk(i-1,j,k,icoff,ijcoff)
135  ijkcm2 = indijk(i-2,j,k,icoff,ijcoff)
136  ijkcp1 = indijk(i+1,j,k,icoff,ijcoff)
137  ijkn = indijk(i ,j,k,inoff,ijnoff)
138  ds = sqrt(si(xcoord,ijkn)*si(xcoord,ijkn)+ &
139  si(ycoord,ijkn)*si(ycoord,ijkn)+ &
140  si(zcoord,ijkn)*si(zcoord,ijkn))
141  nx = si(xcoord,ijkn)/ds
142  ny = si(ycoord,ijkn)/ds
143  nz = si(zcoord,ijkn)/ds
144  svel = sivel(ijkn*indsvel)/ds
145  vola = (0.5_rfreal*(vol(ijkc0)+vol(ijkcm1)))**1.5_rfreal
146  ds = ds*beta5
147 
148 ! ----- limited extrapolations
149 
150  dvar(1) = cv(cv_mixt_dens,ijkc0 ) - cv(cv_mixt_dens,ijkcm1)
151  dvar(2) = dv(dv_mixt_uvel,ijkc0 ) - dv(dv_mixt_uvel,ijkcm1)
152  dvar(3) = dv(dv_mixt_vvel,ijkc0 ) - dv(dv_mixt_vvel,ijkcm1)
153  dvar(4) = dv(dv_mixt_wvel,ijkc0 ) - dv(dv_mixt_wvel,ijkcm1)
154  dvar(5) = dv(dv_mixt_pres,ijkc0 ) - dv(dv_mixt_pres,ijkcm1)
155 
156  dvarm(1) = cv(cv_mixt_dens,ijkcm1) - cv(cv_mixt_dens,ijkcm2)
157  dvarm(2) = dv(dv_mixt_uvel,ijkcm1) - dv(dv_mixt_uvel,ijkcm2)
158  dvarm(3) = dv(dv_mixt_vvel,ijkcm1) - dv(dv_mixt_vvel,ijkcm2)
159  dvarm(4) = dv(dv_mixt_wvel,ijkcm1) - dv(dv_mixt_wvel,ijkcm2)
160  dvarm(5) = dv(dv_mixt_pres,ijkcm1) - dv(dv_mixt_pres,ijkcm2)
161 
162  dvarp(1) = cv(cv_mixt_dens,ijkcp1) - cv(cv_mixt_dens,ijkc0 )
163  dvarp(2) = dv(dv_mixt_uvel,ijkcp1) - dv(dv_mixt_uvel,ijkc0 )
164  dvarp(3) = dv(dv_mixt_vvel,ijkcp1) - dv(dv_mixt_vvel,ijkc0 )
165  dvarp(4) = dv(dv_mixt_wvel,ijkcp1) - dv(dv_mixt_wvel,ijkc0 )
166  dvarp(5) = dv(dv_mixt_pres,ijkcp1) - dv(dv_mixt_pres,ijkc0 )
167 
168  eps2n = eps2(1)*vola
169  deltl(1) = 0.5_rfreal*muscl3(dvar(1) ,dvarm(1),eps2n)
170  deltr(1) = 0.5_rfreal*muscl3(dvarp(1),dvar(1) ,eps2n)
171  eps2n = eps2(2)*vola
172  deltl(2) = 0.5_rfreal*muscl3(dvar(2) ,dvarm(2),eps2n)
173  deltr(2) = 0.5_rfreal*muscl3(dvarp(2),dvar(2) ,eps2n)
174  deltl(3) = 0.5_rfreal*muscl3(dvar(3) ,dvarm(3),eps2n)
175  deltr(3) = 0.5_rfreal*muscl3(dvarp(3),dvar(3) ,eps2n)
176  deltl(4) = 0.5_rfreal*muscl3(dvar(4) ,dvarm(4),eps2n)
177  deltr(4) = 0.5_rfreal*muscl3(dvarp(4),dvar(4) ,eps2n)
178  eps2n = eps2(3)*vola
179  deltl(5) = 0.5_rfreal*muscl3(dvar(5) ,dvarm(5),eps2n)
180  deltr(5) = 0.5_rfreal*muscl3(dvarp(5),dvar(5) ,eps2n)
181 
182 ! ----- left and right states
183 
184  rl = cv(cv_mixt_dens,ijkcm1) + deltl(1)
185  ul = dv(dv_mixt_uvel,ijkcm1) + deltl(2)
186  vl = dv(dv_mixt_vvel,ijkcm1) + deltl(3)
187  wl = dv(dv_mixt_wvel,ijkcm1) + deltl(4)
188  pl = dv(dv_mixt_pres,ijkcm1) + deltl(5)
189  rgas = 8314.3_rfreal/gv(gv_mixt_mol,ijkcm1*indmol)
190  gaml = gv(gv_mixt_cp,ijkcm1*indcp)/(gv(gv_mixt_cp,ijkcm1*indcp)-rgas)
191  ggm1 = gaml/(gaml-1._rfreal)
192  hl = ggm1*pl/rl + 0.5_rfreal*(ul*ul+vl*vl+wl*wl)
193 
194  rr = cv(cv_mixt_dens,ijkc0) - deltr(1)
195  ur = dv(dv_mixt_uvel,ijkc0) - deltr(2)
196  vr = dv(dv_mixt_vvel,ijkc0) - deltr(3)
197  wr = dv(dv_mixt_wvel,ijkc0) - deltr(4)
198  pr = dv(dv_mixt_pres,ijkc0) - deltr(5)
199  rgas = 8314.3_rfreal/gv(gv_mixt_mol,ijkc0*indmol)
200  gamr = gv(gv_mixt_cp,ijkc0*indcp)/(gv(gv_mixt_cp,ijkc0*indcp)-rgas)
201  ggm1 = gamr/(gamr-1._rfreal)
202  hr = ggm1*pr/rr + 0.5_rfreal*(ur*ur+vr*vr+wr*wr)
203 
204 ! ----- Roe`s average
205 
206  rav = sqrt(rl*rr)
207  dd = rav/rl
208  dd1 = 1._rfreal/(1._rfreal+dd)
209  uav = (ul+dd*ur)*dd1
210  vav = (vl+dd*vr)*dd1
211  wav = (wl+dd*wr)*dd1
212  hav = (hl+dd*hr)*dd1
213  q2a = 0.5_rfreal*(uav*uav+vav*vav+wav*wav)
214  gam1 = 0.5_rfreal*(gaml+gamr) - 1._rfreal
215  c2a = gam1*(hav-q2a)
216  cav = sqrt(c2a)
217  uvw = uav*nx + vav*ny + wav*nz - svel
218  du = (ur-ul)*nx + (vr-vl)*ny + (wr-wl)*nz
219 
220 ! ----- eigenvalues
221 
222  h1 = abs(uvw - cav)
223  h2 = abs(uvw)
224  eabs5 = abs(uvw + cav)
225  delta = epsentr*eabs5
226  eabs1 = entropy_corr2( h1,delta )
227  eabs2 = entropy_corr2( h2,delta )
228 
229 ! ----- upwind fluxes
230 
231  h1 = rav*cav*du
232  h2 = eabs1*(pr-pl - h1)/(2._rfreal*c2a)
233  h3 = eabs2*(rr-rl - (pr-pl)/c2a)
234  h4 = eabs2*rav
235  h5 = eabs5*(pr-pl + h1)/(2._rfreal*c2a)
236 
237  fd(1) = h2 + h3 + h5
238  fd(2) = h2*(uav-cav*nx) + h3*uav + h4*(ur-ul-du*nx) + &
239  h5*(uav+cav*nx)
240  fd(3) = h2*(vav-cav*ny) + h3*vav + h4*(vr-vl-du*ny) + &
241  h5*(vav+cav*ny)
242  fd(4) = h2*(wav-cav*nz) + h3*wav + h4*(wr-wl-du*nz) + &
243  h5*(wav+cav*nz)
244  fd(5) = h2*(hav-cav*uvw) + h3*q2a + h4*(uav*(ur-ul)+ &
245  vav*(vr-vl)+wav*(wr-wl)-uvw*du) + &
246  h5*(hav+cav*uvw)
247 
248  fd(1) = fd(1)*ds
249  fd(2) = fd(2)*ds
250  fd(3) = fd(3)*ds
251  fd(4) = fd(4)*ds
252  fd(5) = fd(5)*ds
253 
254  diss(cv_mixt_dens,ijkc0 ) = diss(cv_mixt_dens,ijkc0 ) - fd(1)
255  diss(cv_mixt_xmom,ijkc0 ) = diss(cv_mixt_xmom,ijkc0 ) - fd(2)
256  diss(cv_mixt_ymom,ijkc0 ) = diss(cv_mixt_ymom,ijkc0 ) - fd(3)
257  diss(cv_mixt_zmom,ijkc0 ) = diss(cv_mixt_zmom,ijkc0 ) - fd(4)
258  diss(cv_mixt_ener,ijkc0 ) = diss(cv_mixt_ener,ijkc0 ) - fd(5)
259 
260  diss(cv_mixt_dens,ijkcm1) = diss(cv_mixt_dens,ijkcm1) + fd(1)
261  diss(cv_mixt_xmom,ijkcm1) = diss(cv_mixt_xmom,ijkcm1) + fd(2)
262  diss(cv_mixt_ymom,ijkcm1) = diss(cv_mixt_ymom,ijkcm1) + fd(3)
263  diss(cv_mixt_zmom,ijkcm1) = diss(cv_mixt_zmom,ijkcm1) + fd(4)
264  diss(cv_mixt_ener,ijkcm1) = diss(cv_mixt_ener,ijkcm1) + fd(5)
265  ENDDO ! i
266  ENDDO ! j
267  ENDDO ! k
268 
269 ! dissipation in j-direction --------------------------------------------------
270 
271  DO k=kpcbeg,kpcend
272  DO i=ipcbeg,ipcend
273  DO j=jpcbeg,jpcend+1
274  ijkc0 = indijk(i,j ,k,icoff,ijcoff)
275  ijkcm1 = indijk(i,j-1,k,icoff,ijcoff)
276  ijkcm2 = indijk(i,j-2,k,icoff,ijcoff)
277  ijkcp1 = indijk(i,j+1,k,icoff,ijcoff)
278  ijkn = indijk(i,j ,k,inoff,ijnoff)
279  ds = sqrt(sj(xcoord,ijkn)*sj(xcoord,ijkn)+ &
280  sj(ycoord,ijkn)*sj(ycoord,ijkn)+ &
281  sj(zcoord,ijkn)*sj(zcoord,ijkn))
282  nx = sj(xcoord,ijkn)/ds
283  ny = sj(ycoord,ijkn)/ds
284  nz = sj(zcoord,ijkn)/ds
285  svel = sjvel(ijkn*indsvel)/ds
286  vola = (0.5_rfreal*(vol(ijkc0)+vol(ijkcm1)))**1.5_rfreal
287  ds = ds*beta5
288 
289 ! ----- limited extrapolations
290 
291  dvar(1) = cv(cv_mixt_dens,ijkc0 ) - cv(cv_mixt_dens,ijkcm1)
292  dvar(2) = dv(dv_mixt_uvel,ijkc0 ) - dv(dv_mixt_uvel,ijkcm1)
293  dvar(3) = dv(dv_mixt_vvel,ijkc0 ) - dv(dv_mixt_vvel,ijkcm1)
294  dvar(4) = dv(dv_mixt_wvel,ijkc0 ) - dv(dv_mixt_wvel,ijkcm1)
295  dvar(5) = dv(dv_mixt_pres,ijkc0 ) - dv(dv_mixt_pres,ijkcm1)
296 
297  dvarm(1) = cv(cv_mixt_dens,ijkcm1) - cv(cv_mixt_dens,ijkcm2)
298  dvarm(2) = dv(dv_mixt_uvel,ijkcm1) - dv(dv_mixt_uvel,ijkcm2)
299  dvarm(3) = dv(dv_mixt_vvel,ijkcm1) - dv(dv_mixt_vvel,ijkcm2)
300  dvarm(4) = dv(dv_mixt_wvel,ijkcm1) - dv(dv_mixt_wvel,ijkcm2)
301  dvarm(5) = dv(dv_mixt_pres,ijkcm1) - dv(dv_mixt_pres,ijkcm2)
302 
303  dvarp(1) = cv(cv_mixt_dens,ijkcp1) - cv(cv_mixt_dens,ijkc0 )
304  dvarp(2) = dv(dv_mixt_uvel,ijkcp1) - dv(dv_mixt_uvel,ijkc0 )
305  dvarp(3) = dv(dv_mixt_vvel,ijkcp1) - dv(dv_mixt_vvel,ijkc0 )
306  dvarp(4) = dv(dv_mixt_wvel,ijkcp1) - dv(dv_mixt_wvel,ijkc0 )
307  dvarp(5) = dv(dv_mixt_pres,ijkcp1) - dv(dv_mixt_pres,ijkc0 )
308 
309  eps2n = eps2(1)*vola
310  deltl(1) = 0.5_rfreal*muscl3(dvar(1) ,dvarm(1),eps2n)
311  deltr(1) = 0.5_rfreal*muscl3(dvarp(1),dvar(1) ,eps2n)
312  eps2n = eps2(2)*vola
313  deltl(2) = 0.5_rfreal*muscl3(dvar(2) ,dvarm(2),eps2n)
314  deltr(2) = 0.5_rfreal*muscl3(dvarp(2),dvar(2) ,eps2n)
315  deltl(3) = 0.5_rfreal*muscl3(dvar(3) ,dvarm(3),eps2n)
316  deltr(3) = 0.5_rfreal*muscl3(dvarp(3),dvar(3) ,eps2n)
317  deltl(4) = 0.5_rfreal*muscl3(dvar(4) ,dvarm(4),eps2n)
318  deltr(4) = 0.5_rfreal*muscl3(dvarp(4),dvar(4) ,eps2n)
319  eps2n = eps2(3)*vola
320  deltl(5) = 0.5_rfreal*muscl3(dvar(5) ,dvarm(5),eps2n)
321  deltr(5) = 0.5_rfreal*muscl3(dvarp(5),dvar(5) ,eps2n)
322 
323 ! ----- left and right states
324 
325  rl = cv(cv_mixt_dens,ijkcm1) + deltl(1)
326  ul = dv(dv_mixt_uvel,ijkcm1) + deltl(2)
327  vl = dv(dv_mixt_vvel,ijkcm1) + deltl(3)
328  wl = dv(dv_mixt_wvel,ijkcm1) + deltl(4)
329  pl = dv(dv_mixt_pres,ijkcm1) + deltl(5)
330  rgas = 8314.3_rfreal/gv(gv_mixt_mol,ijkcm1*indmol)
331  gaml = gv(gv_mixt_cp,ijkcm1*indcp)/(gv(gv_mixt_cp,ijkcm1*indcp)-rgas)
332  ggm1 = gaml/(gaml-1._rfreal)
333  hl = ggm1*pl/rl + 0.5_rfreal*(ul*ul+vl*vl+wl*wl)
334 
335  rr = cv(cv_mixt_dens,ijkc0) - deltr(1)
336  ur = dv(dv_mixt_uvel,ijkc0) - deltr(2)
337  vr = dv(dv_mixt_vvel,ijkc0) - deltr(3)
338  wr = dv(dv_mixt_wvel,ijkc0) - deltr(4)
339  pr = dv(dv_mixt_pres,ijkc0) - deltr(5)
340  rgas = 8314.3_rfreal/gv(gv_mixt_mol,ijkc0*indmol)
341  gamr = gv(gv_mixt_cp,ijkc0*indcp)/(gv(gv_mixt_cp,ijkc0*indcp)-rgas)
342  ggm1 = gamr/(gamr-1._rfreal)
343  hr = ggm1*pr/rr + 0.5_rfreal*(ur*ur+vr*vr+wr*wr)
344 
345 ! ----- Roe`s average
346 
347  rav = sqrt(rl*rr)
348  dd = rav/rl
349  dd1 = 1._rfreal/(1._rfreal+dd)
350  uav = (ul+dd*ur)*dd1
351  vav = (vl+dd*vr)*dd1
352  wav = (wl+dd*wr)*dd1
353  hav = (hl+dd*hr)*dd1
354  q2a = 0.5_rfreal*(uav*uav+vav*vav+wav*wav)
355  gam1 = 0.5_rfreal*(gaml+gamr) - 1._rfreal
356  c2a = gam1*(hav-q2a)
357  cav = sqrt(c2a)
358  uvw = uav*nx + vav*ny + wav*nz - svel
359  du = (ur-ul)*nx + (vr-vl)*ny + (wr-wl)*nz
360 
361 ! ----- eigenvalues
362 
363  h1 = abs(uvw - cav)
364  h2 = abs(uvw)
365  eabs5 = abs(uvw + cav)
366  delta = epsentr*eabs5
367  eabs1 = entropy_corr2( h1,delta )
368  eabs2 = entropy_corr2( h2,delta )
369 
370 ! ----- upwind fluxes
371 
372  h1 = rav*cav*du
373  h2 = eabs1*(pr-pl - h1)/(2._rfreal*c2a)
374  h3 = eabs2*(rr-rl - (pr-pl)/c2a)
375  h4 = eabs2*rav
376  h5 = eabs5*(pr-pl + h1)/(2._rfreal*c2a)
377 
378  fd(1) = h2 + h3 + h5
379  fd(2) = h2*(uav-cav*nx) + h3*uav + h4*(ur-ul-du*nx) + &
380  h5*(uav+cav*nx)
381  fd(3) = h2*(vav-cav*ny) + h3*vav + h4*(vr-vl-du*ny) + &
382  h5*(vav+cav*ny)
383  fd(4) = h2*(wav-cav*nz) + h3*wav + h4*(wr-wl-du*nz) + &
384  h5*(wav+cav*nz)
385  fd(5) = h2*(hav-cav*uvw) + h3*q2a + h4*(uav*(ur-ul)+ &
386  vav*(vr-vl)+wav*(wr-wl)-uvw*du) + &
387  h5*(hav+cav*uvw)
388 
389  fd(1) = fd(1)*ds
390  fd(2) = fd(2)*ds
391  fd(3) = fd(3)*ds
392  fd(4) = fd(4)*ds
393  fd(5) = fd(5)*ds
394 
395  diss(cv_mixt_dens,ijkc0 ) = diss(cv_mixt_dens,ijkc0 ) - fd(1)
396  diss(cv_mixt_xmom,ijkc0 ) = diss(cv_mixt_xmom,ijkc0 ) - fd(2)
397  diss(cv_mixt_ymom,ijkc0 ) = diss(cv_mixt_ymom,ijkc0 ) - fd(3)
398  diss(cv_mixt_zmom,ijkc0 ) = diss(cv_mixt_zmom,ijkc0 ) - fd(4)
399  diss(cv_mixt_ener,ijkc0 ) = diss(cv_mixt_ener,ijkc0 ) - fd(5)
400 
401  diss(cv_mixt_dens,ijkcm1) = diss(cv_mixt_dens,ijkcm1) + fd(1)
402  diss(cv_mixt_xmom,ijkcm1) = diss(cv_mixt_xmom,ijkcm1) + fd(2)
403  diss(cv_mixt_ymom,ijkcm1) = diss(cv_mixt_ymom,ijkcm1) + fd(3)
404  diss(cv_mixt_zmom,ijkcm1) = diss(cv_mixt_zmom,ijkcm1) + fd(4)
405  diss(cv_mixt_ener,ijkcm1) = diss(cv_mixt_ener,ijkcm1) + fd(5)
406  ENDDO ! j
407  ENDDO ! i
408  ENDDO ! k
409 
410 ! dissipation in k-direction --------------------------------------------------
411 
412  DO j=jpcbeg,jpcend
413  DO i=ipcbeg,ipcend
414  DO k=kpcbeg,kpcend+1
415  ijkc0 = indijk(i,j,k ,icoff,ijcoff)
416  ijkcm1 = indijk(i,j,k-1,icoff,ijcoff)
417  ijkcm2 = indijk(i,j,k-2,icoff,ijcoff)
418  ijkcp1 = indijk(i,j,k+1,icoff,ijcoff)
419  ijkn = indijk(i,j,k ,inoff,ijnoff)
420  ds = sqrt(sk(xcoord,ijkn)*sk(xcoord,ijkn)+ &
421  sk(ycoord,ijkn)*sk(ycoord,ijkn)+ &
422  sk(zcoord,ijkn)*sk(zcoord,ijkn))
423  nx = sk(xcoord,ijkn)/ds
424  ny = sk(ycoord,ijkn)/ds
425  nz = sk(zcoord,ijkn)/ds
426  svel = skvel(ijkn*indsvel)/ds
427  vola = (0.5_rfreal*(vol(ijkc0)+vol(ijkcm1)))**1.5_rfreal
428  ds = ds*beta5
429 
430 ! ----- limited extrapolations
431 
432  dvar(1) = cv(cv_mixt_dens,ijkc0 ) - cv(cv_mixt_dens,ijkcm1)
433  dvar(2) = dv(dv_mixt_uvel,ijkc0 ) - dv(dv_mixt_uvel,ijkcm1)
434  dvar(3) = dv(dv_mixt_vvel,ijkc0 ) - dv(dv_mixt_vvel,ijkcm1)
435  dvar(4) = dv(dv_mixt_wvel,ijkc0 ) - dv(dv_mixt_wvel,ijkcm1)
436  dvar(5) = dv(dv_mixt_pres,ijkc0 ) - dv(dv_mixt_pres,ijkcm1)
437 
438  dvarm(1) = cv(cv_mixt_dens,ijkcm1) - cv(cv_mixt_dens,ijkcm2)
439  dvarm(2) = dv(dv_mixt_uvel,ijkcm1) - dv(dv_mixt_uvel,ijkcm2)
440  dvarm(3) = dv(dv_mixt_vvel,ijkcm1) - dv(dv_mixt_vvel,ijkcm2)
441  dvarm(4) = dv(dv_mixt_wvel,ijkcm1) - dv(dv_mixt_wvel,ijkcm2)
442  dvarm(5) = dv(dv_mixt_pres,ijkcm1) - dv(dv_mixt_pres,ijkcm2)
443 
444  dvarp(1) = cv(cv_mixt_dens,ijkcp1) - cv(cv_mixt_dens,ijkc0 )
445  dvarp(2) = dv(dv_mixt_uvel,ijkcp1) - dv(dv_mixt_uvel,ijkc0 )
446  dvarp(3) = dv(dv_mixt_vvel,ijkcp1) - dv(dv_mixt_vvel,ijkc0 )
447  dvarp(4) = dv(dv_mixt_wvel,ijkcp1) - dv(dv_mixt_wvel,ijkc0 )
448  dvarp(5) = dv(dv_mixt_pres,ijkcp1) - dv(dv_mixt_pres,ijkc0 )
449 
450  eps2n = eps2(1)*vola
451  deltl(1) = 0.5_rfreal*muscl3(dvar(1) ,dvarm(1),eps2n)
452  deltr(1) = 0.5_rfreal*muscl3(dvarp(1),dvar(1) ,eps2n)
453  eps2n = eps2(2)*vola
454  deltl(2) = 0.5_rfreal*muscl3(dvar(2) ,dvarm(2),eps2n)
455  deltr(2) = 0.5_rfreal*muscl3(dvarp(2),dvar(2) ,eps2n)
456  deltl(3) = 0.5_rfreal*muscl3(dvar(3) ,dvarm(3),eps2n)
457  deltr(3) = 0.5_rfreal*muscl3(dvarp(3),dvar(3) ,eps2n)
458  deltl(4) = 0.5_rfreal*muscl3(dvar(4) ,dvarm(4),eps2n)
459  deltr(4) = 0.5_rfreal*muscl3(dvarp(4),dvar(4) ,eps2n)
460  eps2n = eps2(3)*vola
461  deltl(5) = 0.5_rfreal*muscl3(dvar(5) ,dvarm(5),eps2n)
462  deltr(5) = 0.5_rfreal*muscl3(dvarp(5),dvar(5) ,eps2n)
463 
464 ! ----- left and right states
465 
466  rl = cv(cv_mixt_dens,ijkcm1) + deltl(1)
467  ul = dv(dv_mixt_uvel,ijkcm1) + deltl(2)
468  vl = dv(dv_mixt_vvel,ijkcm1) + deltl(3)
469  wl = dv(dv_mixt_wvel,ijkcm1) + deltl(4)
470  pl = dv(dv_mixt_pres,ijkcm1) + deltl(5)
471  rgas = 8314.3_rfreal/gv(gv_mixt_mol,ijkcm1*indmol)
472  gaml = gv(gv_mixt_cp,ijkcm1*indcp)/(gv(gv_mixt_cp,ijkcm1*indcp)-rgas)
473  ggm1 = gaml/(gaml-1._rfreal)
474  hl = ggm1*pl/rl + 0.5_rfreal*(ul*ul+vl*vl+wl*wl)
475 
476  rr = cv(cv_mixt_dens,ijkc0) - deltr(1)
477  ur = dv(dv_mixt_uvel,ijkc0) - deltr(2)
478  vr = dv(dv_mixt_vvel,ijkc0) - deltr(3)
479  wr = dv(dv_mixt_wvel,ijkc0) - deltr(4)
480  pr = dv(dv_mixt_pres,ijkc0) - deltr(5)
481  rgas = 8314.3_rfreal/gv(gv_mixt_mol,ijkc0*indmol)
482  gamr = gv(gv_mixt_cp,ijkc0*indcp)/(gv(gv_mixt_cp,ijkc0*indcp)-rgas)
483  ggm1 = gamr/(gamr-1._rfreal)
484  hr = ggm1*pr/rr + 0.5_rfreal*(ur*ur+vr*vr+wr*wr)
485 
486 ! ----- Roe`s average
487 
488  rav = sqrt(rl*rr)
489  dd = rav/rl
490  dd1 = 1._rfreal/(1._rfreal+dd)
491  uav = (ul+dd*ur)*dd1
492  vav = (vl+dd*vr)*dd1
493  wav = (wl+dd*wr)*dd1
494  hav = (hl+dd*hr)*dd1
495  q2a = 0.5_rfreal*(uav*uav+vav*vav+wav*wav)
496  gam1 = 0.5_rfreal*(gaml+gamr) - 1._rfreal
497  c2a = gam1*(hav-q2a)
498  cav = sqrt(c2a)
499  uvw = uav*nx + vav*ny + wav*nz - svel
500  du = (ur-ul)*nx + (vr-vl)*ny + (wr-wl)*nz
501 
502 ! ----- eigenvalues
503 
504  h1 = abs(uvw - cav)
505  h2 = abs(uvw)
506  eabs5 = abs(uvw + cav)
507  delta = epsentr*eabs5
508  eabs1 = entropy_corr2( h1,delta )
509  eabs2 = entropy_corr2( h2,delta )
510 
511 ! ----- upwind fluxes
512 
513  h1 = rav*cav*du
514  h2 = eabs1*(pr-pl - h1)/(2._rfreal*c2a)
515  h3 = eabs2*(rr-rl - (pr-pl)/c2a)
516  h4 = eabs2*rav
517  h5 = eabs5*(pr-pl + h1)/(2._rfreal*c2a)
518 
519  fd(1) = h2 + h3 + h5
520  fd(2) = h2*(uav-cav*nx) + h3*uav + h4*(ur-ul-du*nx) + &
521  h5*(uav+cav*nx)
522  fd(3) = h2*(vav-cav*ny) + h3*vav + h4*(vr-vl-du*ny) + &
523  h5*(vav+cav*ny)
524  fd(4) = h2*(wav-cav*nz) + h3*wav + h4*(wr-wl-du*nz) + &
525  h5*(wav+cav*nz)
526  fd(5) = h2*(hav-cav*uvw) + h3*q2a + h4*(uav*(ur-ul)+ &
527  vav*(vr-vl)+wav*(wr-wl)-uvw*du) + &
528  h5*(hav+cav*uvw)
529 
530  fd(1) = fd(1)*ds
531  fd(2) = fd(2)*ds
532  fd(3) = fd(3)*ds
533  fd(4) = fd(4)*ds
534  fd(5) = fd(5)*ds
535 
536  diss(cv_mixt_dens,ijkc0 ) = diss(cv_mixt_dens,ijkc0 ) - fd(1)
537  diss(cv_mixt_xmom,ijkc0 ) = diss(cv_mixt_xmom,ijkc0 ) - fd(2)
538  diss(cv_mixt_ymom,ijkc0 ) = diss(cv_mixt_ymom,ijkc0 ) - fd(3)
539  diss(cv_mixt_zmom,ijkc0 ) = diss(cv_mixt_zmom,ijkc0 ) - fd(4)
540  diss(cv_mixt_ener,ijkc0 ) = diss(cv_mixt_ener,ijkc0 ) - fd(5)
541 
542  diss(cv_mixt_dens,ijkcm1) = diss(cv_mixt_dens,ijkcm1) + fd(1)
543  diss(cv_mixt_xmom,ijkcm1) = diss(cv_mixt_xmom,ijkcm1) + fd(2)
544  diss(cv_mixt_ymom,ijkcm1) = diss(cv_mixt_ymom,ijkcm1) + fd(3)
545  diss(cv_mixt_zmom,ijkcm1) = diss(cv_mixt_zmom,ijkcm1) + fd(4)
546  diss(cv_mixt_ener,ijkcm1) = diss(cv_mixt_ener,ijkcm1) + fd(5)
547  ENDDO ! k
548  ENDDO ! i
549  ENDDO ! j
550 
551 ! finalize --------------------------------------------------------------------
552 
553  CALL deregisterfunction( region%global )
554 
555 ! #############################################################################
556 
557  CONTAINS
558 
559  REAL(RFREAL) FUNCTION entropy_corr2( z,d )
560 
561  REAL(RFREAL) :: z, d
562 
563  IF (z > d) THEN
564  entropy_corr2 = z
565  ELSE
566  entropy_corr2 = 0.5_rfreal*(z*z+d*d)/d
567  ENDIF
568 
569  END FUNCTION entropy_corr2
570 
571 END SUBROUTINE rflo_roedissipsecond
572 
573 !******************************************************************************
574 !
575 ! RCS Revision history:
576 !
577 ! $Log: RFLO_RoeDissipSecond.F90,v $
578 ! Revision 1.3 2008/12/06 08:44:27 mtcampbe
579 ! Updated license.
580 !
581 ! Revision 1.2 2008/11/19 22:17:38 mtcampbe
582 ! Added Illinois Open Source License/Copyright
583 !
584 ! Revision 1.1 2004/11/29 20:51:40 wasistho
585 ! lower to upper case
586 !
587 ! Revision 1.6 2003/11/20 16:40:40 mdbrandy
588 ! Backing out RocfluidMP changes from 11-17-03
589 !
590 ! Revision 1.2 2003/10/01 23:52:10 jblazek
591 ! Corrected bug in moving noslip wall BC and grid speeds.
592 !
593 ! Revision 1.1 2003/05/29 17:28:43 jblazek
594 ! Implemented Roe scheme.
595 !
596 !******************************************************************************
597 
598 
599 
600 
601 
602 
603 
d
const NT & d
Definition: rational_rotation.h:73
k
j indices k indices k
Definition: Indexing.h:6
kpcbeg
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE kpcbeg
Definition: RFLO_GetDimensPhys.h:27
moderror::registerfunction
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
rflo_roedissipsecond
subroutine rflo_roedissipsecond(region)
Definition: RFLO_RoeDissipSecond.F90:44
sqrt
double sqrt(double d)
Definition: double.h:73
jpcbeg
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE jpcbeg
Definition: RFLO_GetDimensPhys.h:27
ipcend
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE ipcend
Definition: RFLO_GetDimensPhys.h:27
rflo_getnodeoffset
subroutine rflo_getnodeoffset(region, iLev, iNodeOffset, ijNodeOffset)
Definition: RFLO_GetDimensions.F90:267
z
void int int int REAL REAL REAL * z
Definition: write.cpp:76
ipcbeg
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE ipcbeg
Definition: RFLO_GetDimensPhys.h:27
i
blockLoc i
Definition: read.cpp:79
rflo_getcelloffset
subroutine rflo_getcelloffset(region, iLev, iCellOffset, ijCellOffset)
Definition: RFLO_GetDimensions.F90:235
Direction_2::delta
RT delta(int i) const
Definition: Direction_2.h:159
j
j indices j
Definition: Indexing.h:6
jpcend
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE jpcend
Definition: RFLO_GetDimensPhys.h:27
moderror::deregisterfunction
subroutine deregisterfunction(global)
Definition: ModError.F90:469
rflo_getdimensphys
subroutine rflo_getdimensphys(region, iLev, ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend)
Definition: RFLO_GetDimensions.F90:109
entropy_corr2
REAL(RFREAL) function entropy_corr2(z, d)
Definition: RFLO_RoeDissipSecond.F90:559