Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLU_CentralSecondPatch_GL.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 central convective fluxes for mixture using second-order
26 ! accurate approximation through a patch by using an average of variables.
27 !
28 ! Description: None.
29 !
30 ! Input:
31 ! pRegion Pointer to region
32 ! pPatch Pointer to patch
33 !
34 ! Output: None.
35 !
36 ! Notes: None.
37 !
38 ! ******************************************************************************
39 !
40 ! $Id: RFLU_CentralSecondPatch_GL.F90,v 1.5 2008/12/06 08:44:29 mtcampbe Exp $
41 !
42 ! Copyright: (c) 2006 by the University of Illinois
43 !
44 ! ******************************************************************************
45 
46 SUBROUTINE rflu_centralsecondpatch_gl(pRegion,pPatch)
47 
48  USE modglobal, ONLY: t_global
49  USE moddatatypes
50  USE modbndpatch, ONLY: t_patch
51  USE moddatastruct, ONLY: t_region
52  USE moderror
53  USE modparameters
54  USE modmpi
55 
60 
61  USE modinterfaces, ONLY: bcondinflowperf_gl, &
63  mixtgasliq_c, &
64  mixtgasliq_p, &
65  mixtliq_c2_bp, &
69  mixtperf_d_prt, &
73  mixtperf_p_drt, &
74  mixtperf_r_m, &
78 
79  IMPLICIT NONE
80 
81 ! ******************************************************************************
82 ! Definitions and declarations
83 ! ******************************************************************************
84 
85 ! ==============================================================================
86 ! Arguments
87 ! ==============================================================================
88 
89  TYPE(t_patch), POINTER :: ppatch
90  TYPE(t_region), POINTER :: pregion
91 
92 ! ==============================================================================
93 ! Locals
94 ! ==============================================================================
95 
96  LOGICAL :: corrflag,decideprintflag
97  CHARACTER(CHRLEN) :: rcsidentstring
98  INTEGER, PARAMETER :: max_inoutflow_locs = 10
99  INTEGER :: c1,bcoptfixed,bcopttype,distrib,ifl,indcp,indgs,indmf,indsd, &
100  indmol,nlocs
101  INTEGER :: loc(max_inoutflow_locs,min_val:max_val)
102  REAL(RFREAL) :: bgl2,bll2,bp,bt,bvl2,cgl2,cll2,cml,cp,cvg,cvl,cvl2, &
103  cvm,cvv,dx,dy,dz,fs,fsu,icpref,irl,mf,mm,nm,nx,ny,nz, &
104  pl,po,pr,pref,press, ql,qr,rl,rel,rer,rg,rgl,rgpgl, &
105  rgpgr,rhogr,rholr,rhovr,rll,rlpll,ro,rr,rref,rul,rur, &
106  rv,rvl,rvpvl,rvpvr,rvr,rwl,rwr, rygl,ryll,ryvl,temp,tl, &
107  to,tr,ul,ur,vfgr,vflr,vfvr,vl,vl2,vm2,vr,vref,wl,wr,xc, &
108  yc,zc,ygl,ygr,yll,ylr,yvl,yvr,vfgl,vfll,vfvl
109  REAL(RFREAL) :: flx(5)
110  REAL(RFREAL), DIMENSION(:), POINTER :: mfmixt
111  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv,dv,gv,pcvspec,rhs,vals,sd
112  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad,pgradspec
113  TYPE(t_global), POINTER :: global
114 
115 ! ******************************************************************************
116 ! Start
117 ! ******************************************************************************
118 
119  rcsidentstring = '$RCSfile: RFLU_CentralSecondPatch_GL.F90,v $'
120 
121  global => pregion%global
122 
123  CALL registerfunction(global,'RFLU_CentralSecondPatch_GL',&
124  'RFLU_CentralSecondPatch_GL.F90')
125 
126 ! ******************************************************************************
127 ! Set pointers and variables
128 ! ******************************************************************************
129 
130  indcp = pregion%mixtInput%indCp
131  indmf = pregion%mixtInput%indMfMixt
132  indsd = pregion%mixtInput%indSd
133  indmol = pregion%mixtInput%indMol
134 
135  indgs = pregion%grid%indGs
136 
137  distrib = ppatch%mixt%distrib
138 
139  cv => pregion%mixt%cv
140  dv => pregion%mixt%dv
141  gv => pregion%mixt%gv
142  grad => pregion%mixt%gradCell
143  rhs => pregion%mixt%rhs
144  sd => pregion%mixt%sd
145 
146  pcvspec => pregion%spec%cv
147  pgradspec => pregion%spec%gradCell
148 
149  mfmixt => ppatch%mfMixt
150 
151  nlocs = 0
152 
153  decideprintflag = rflu_decideprint(global)
154 
155  pref = global%refPressure
156  rref = global%refDensity
157  vref = global%refVelocity
158 
159  icpref = 2.0_rfreal/(rref*vref*vref)
160 
161 ! ******************************************************************************
162 ! Define constants
163 ! ******************************************************************************
164 
165  ro = global%refDensityLiq
166  po = global%refPressLiq
167  to = global%refTempLiq
168  bp = global%refBetaPLiq
169  bt = global%refBetaTLiq
170  cvl = global%refCvLiq
171 
172  rg = mixtperf_r_m(pregion%specInput%specType(1)%pMaterial%molw)
173  cvg = mixtperf_cv_cpr(pregion%specInput%specType(1)%pMaterial%spht,rg)
174 
175  rv = mixtperf_r_m(pregion%specInput%specType(2)%pMaterial%molw)
176  cvv = mixtperf_cv_cpr(pregion%specInput%specType(2)%pMaterial%spht,rv)
177 
178 ! ******************************************************************************
179 ! Select boundary type
180 ! ******************************************************************************
181 
182  SELECT CASE ( ppatch%bcType )
183 
184 ! ==============================================================================
185 ! Inflow based on velocity and temperature
186 ! ==============================================================================
187 
188  CASE ( bc_inflow_veltemp )
189  vals => ppatch%mixt%vals
190  bcopttype = ppatch%mixt%switches(bcswi_inflow_type)
191 
192  DO ifl = 1,ppatch%nBFaces
193  c1 = ppatch%bf2c(ifl)
194 
195  nx = ppatch%fn(xcoord,ifl)
196  ny = ppatch%fn(ycoord,ifl)
197  nz = ppatch%fn(zcoord,ifl)
198  nm = ppatch%fn(xyzmag,ifl)
199 
200  xc = ppatch%fc(xcoord,ifl)
201  yc = ppatch%fc(ycoord,ifl)
202  zc = ppatch%fc(zcoord,ifl)
203 
204  fs = ppatch%gs(indgs*ifl)
205 
206  rl = cv(cv_mixt_dens,c1)
207  irl = 1.0_rfreal/rl
208 
209  ul = cv(cv_mixt_xmom,c1)*irl
210  vl = cv(cv_mixt_ymom,c1)*irl
211  wl = cv(cv_mixt_zmom,c1)*irl
212  tl = dv(dv_mixt_temp,c1)
213 
214  ygl = pcvspec(1,c1)*irl
215  yvl = pcvspec(2,c1)*irl
216 
217  dx = xc - pregion%grid%cofg(xcoord,c1)
218  dy = yc - pregion%grid%cofg(ycoord,c1)
219  dz = zc - pregion%grid%cofg(zcoord,c1)
220 
221  rl = rl + grad(xcoord,grc_mixt_dens,c1)*dx &
222  + grad(ycoord,grc_mixt_dens,c1)*dy &
223  + grad(zcoord,grc_mixt_dens,c1)*dz
224  ul = ul + grad(xcoord,grc_mixt_xvel,c1)*dx &
225  + grad(ycoord,grc_mixt_xvel,c1)*dy &
226  + grad(zcoord,grc_mixt_xvel,c1)*dz
227  vl = vl + grad(xcoord,grc_mixt_yvel,c1)*dx &
228  + grad(ycoord,grc_mixt_yvel,c1)*dy &
229  + grad(zcoord,grc_mixt_yvel,c1)*dz
230  wl = wl + grad(xcoord,grc_mixt_zvel,c1)*dx &
231  + grad(ycoord,grc_mixt_zvel,c1)*dy &
232  + grad(zcoord,grc_mixt_zvel,c1)*dz
233  tl = tl + grad(xcoord,grc_mixt_temp,c1)*dx &
234  + grad(ycoord,grc_mixt_temp,c1)*dy &
235  + grad(zcoord,grc_mixt_temp,c1)*dz
236 
237  ygl = ygl + pgradspec(xcoord,1,c1)*dx &
238  + pgradspec(ycoord,1,c1)*dy &
239  + pgradspec(zcoord,1,c1)*dz
240  yvl = yvl + pgradspec(xcoord,2,c1)*dx &
241  + pgradspec(ycoord,2,c1)*dy &
242  + pgradspec(zcoord,2,c1)*dz
243 
244  IF ( ygl > 1.0_rfreal) THEN
245  ygl = 1.0_rfreal
246  yvl = 0.0_rfreal
247  ELSE IF ( ygl < 0.0_rfreal ) THEN
248  ygl = 0.0_rfreal
249  ELSE IF ( yvl > 1.0_rfreal ) THEN
250  yvl = 1.0_rfreal
251  ygl = 0.0_rfreal
252  ELSE IF ( yvl < 0.0_rfreal ) THEN
253  yvl = 0.0_rfreal
254  END IF ! Ygl
255 
256  yll = 1.0_rfreal - ygl - yvl
257 
258  IF ( yll > 1.0_rfreal) THEN
259  yll = 1.0_rfreal
260  ELSE IF ( yll < 0.0_rfreal) THEN
261  yll = 0.0_rfreal
262  END IF ! Yll
263 
264  vl2 = ul*ul + vl*vl + wl*wl
265 
266  cll2 = mixtliq_c2_bp(bp)
267  cvl2 = mixtperf_c2_grt(1.0_rfreal,rv,tl)
268  cgl2 = mixtperf_c2_grt(1.0_rfreal,rg,tl)
269 
270  rygl = rl*ygl
271  ryvl = rl*yvl
272  ryll = rl*yll
273 
274  pl = mixtgasliq_p(ryll,ryvl,rygl,cll2,cvl2,cgl2,rl, &
275  ro,po,to,bp,bt,tl)
276 
277  rll = mixtliq_d_dobpppobttto(ro,bp,bt,pl,po,tl,to)
278  rvl = mixtperf_d_prt(pl,rv,tl)
279  rgl = mixtperf_d_prt(pl,rg,tl)
280 
281  vfgl = rygl/rgl
282  vfvl = ryvl/rvl
283  vfll = ryll/rll
284 
285  cvm = (ryll*cvl + ryvl*cvv + rygl*cvg)/rl
286 
287  bll2 = -bt/bp
288  bvl2 = rvl*rv
289  bgl2 = rgl*rg
290 
291  cml = mixtgasliq_c(cvm,rl,pl,rll,rvl,rgl,vfll,vfvl,vfgl, &
292  cll2,cvl2,cgl2,bll2,bvl2,bgl2)
293  rul = rl*ul
294  rvl = rl*vl
295  rwl = rl*wl
296  rel = rl*cvm*tl + 0.5*rl*vl2
297 
298  rgpgl = rygl
299  rvpvl = ryvl
300  rlpll = ryll
301 
302  ur = vals(bcdat_inflow_u,distrib*ifl)
303  vr = vals(bcdat_inflow_v,distrib*ifl)
304  wr = vals(bcdat_inflow_w,distrib*ifl)
305  temp = vals(bcdat_inflow_t,distrib*ifl)
306 
307  ygr = ppatch%spec%vals(1,distrib*ifl)
308  yvr = ppatch%spec%vals(2,distrib*ifl)
309  ylr = 1.0_rfreal - ygr - yvr
310 
311  IF ( bcopttype /= bcopt_subsonic ) THEN
312  press = vals(bcdat_inflow_p, distrib*ifl)
313  rholr = mixtliq_d_dobpppobttto(ro,bp,bt,press,po,temp,to)
314  rhogr = mixtperf_d_prt(press,rg,temp)
315  rhovr = mixtperf_d_prt(press,rv,temp)
316  rr = 1.0_rfreal/(ygr/rhogr + yvr/rhovr + ylr/rholr)
317  vfgr = (rr*ygr)/rhogr
318  vfvr = (rr*yvr)/rhovr
319  vflr = (rr*ylr)/rholr
320  ELSE
321  press = pl+(sqrt(ur*ur+vr*vr+wr*wr) &
322  - sqrt(ul*ul+vl*vl+wl*wl))*rl*cml
323  rholr = mixtliq_d_dobpppobttto(ro,bp,bt,press,po,temp,to)
324  rhogr = mixtperf_d_prt(press,rg,temp)
325  rhovr = mixtperf_d_prt(press,rv,temp)
326  rr = 1.0_rfreal/(ygr/rhogr + yvr/rhovr + ylr/rholr)
327  vfgr = (rr*ygr)/rhogr
328  vfvr = (rr*yvr)/rhovr
329  vflr = (rr*ylr)/rholr
330  END IF ! bcOptType
331 
332  CALL bcondinflowperf_gl(bcopttype,ro,po,to,bp,bt,cvl,cvv,cvg,rg,rv, &
333  ur,vr,wr,vfgr,vfvr,vflr,temp,press,nx,ny, &
334  nz,rl,rul,rvl,rwl,rel,rgpgl,rvpvl,pl,rr, &
335  rur,rvr,rwr,rer,rgpgr,rvpvr,pr)
336 
337  ql = (rul*nx + rvl*ny + rwl*nz)/rl - fs
338  qr = (rur*nx + rvr*ny + rwr*nz)/rr - fs
339 
340  flx(1) = 0.5_rfreal*(ql* rl + qr* rr )*nm
341  flx(2) = 0.5_rfreal*(ql* rul + pl*nx + qr* rur + pr*nx)*nm
342  flx(3) = 0.5_rfreal*(ql* rvl + pl*ny + qr* rvr + pr*ny)*nm
343  flx(4) = 0.5_rfreal*(ql* rwl + pl*nz + qr* rwr + pr*nz)*nm
344  flx(5) = 0.5_rfreal*(ql*(rel + pl) + pl*fs + qr*(rer + pr) + pr*fs)*nm
345 
346  mfmixt(indmf*ifl) = flx(1)
347  ppatch%cp(ifl) = icpref*(0.5_rfreal*(pl + pr) - pref)
348 
349  rhs(cv_mixt_dens,c1) = rhs(cv_mixt_dens,c1) + flx(1)
350  rhs(cv_mixt_xmom,c1) = rhs(cv_mixt_xmom,c1) + flx(2)
351  rhs(cv_mixt_ymom,c1) = rhs(cv_mixt_ymom,c1) + flx(3)
352  rhs(cv_mixt_zmom,c1) = rhs(cv_mixt_zmom,c1) + flx(4)
353  rhs(cv_mixt_ener,c1) = rhs(cv_mixt_ener,c1) + flx(5)
354 
355  sd(sd_xmom,c1*indsd) = sd(sd_xmom,c1*indsd) &
356  + 0.5_rfreal*(rul/rl + rur/rr)*flx(1)
357  sd(sd_ymom,c1*indsd) = sd(sd_ymom,c1*indsd) &
358  + 0.5_rfreal*(rvl/rl + rvr/rr)*flx(1)
359  sd(sd_zmom,c1*indsd) = sd(sd_zmom,c1*indsd) &
360  + 0.5_rfreal*(rwl/rl + rwr/rr)*flx(1)
361 
362  IF ( pregion%irkStep == 1 ) THEN
363  global%massIn = global%massIn - min(flx(1),0.0_rfreal)
364  global%massOut = global%massOut + max(flx(1),0.0_rfreal)
365  END IF ! pRegion
366 
367  IF ( (global%checkLevel == check_high) .AND. &
368  (global%verbLevel >= verbose_high) .AND. &
369  (global%myProcid == masterproc) .AND. &
370  (decideprintflag .EQV. .true.) ) THEN
371  IF ( flx(1) > 0.0_rfreal ) THEN
372  nlocs = nlocs + 1
373 
374  IF ( nlocs == 1 ) THEN
375  global%warnCounter = global%warnCounter + 1
376 
377  WRITE(stdout,'(A,1X,A,1X,I9)') solver_name, &
378  '*** WARNING *** Outflow detected at inflow boundary!'
379  WRITE(stdout,'(A,3X,A,1X,I5.5)') solver_name,'Global region:', &
380  pregion%iRegionGlobal
381  IF ( global%flowType == flow_unsteady ) THEN
382  WRITE(stdout,'(A,3X,A,1X,1PE11.5)') solver_name, &
383  'Current time:', &
384  global%currentTime
385  ELSE
386  WRITE(stdout,'(A,3X,A,1X,I6.6)') solver_name, &
387  'Current iteration number:', &
388  global%currentIter
389  END IF ! global%flowType
390 
391  WRITE(stdout,'(A,3X,A,1X,I1)') solver_name, &
392  'Runge-Kutta stage:', &
393  pregion%irkStep
394  END IF ! nLocs
395 
396  IF ( nlocs <= max_inoutflow_locs ) THEN
397  loc(nlocs,min_val:max_val) = c1
398  END IF ! nLocs
399  END IF ! flx(1)
400  END IF ! global%checkLevel
401  END DO ! ifl
402 
403 ! ==============================================================================
404 ! Outflow
405 ! ==============================================================================
406 
407  CASE ( bc_outflow )
408  bcopttype = ppatch%mixt%switches(bcswi_outflow_type)
409 
410  IF ( bcopttype /= bcopt_supersonic ) THEN
411  vals => ppatch%mixt%vals
412  END IF ! bcOptType
413 
414  DO ifl = 1,ppatch%nBFaces
415  c1 = ppatch%bf2c(ifl)
416 
417  nx = ppatch%fn(xcoord,ifl)
418  ny = ppatch%fn(ycoord,ifl)
419  nz = ppatch%fn(zcoord,ifl)
420  nm = ppatch%fn(xyzmag,ifl)
421 
422  xc = ppatch%fc(xcoord,ifl)
423  yc = ppatch%fc(ycoord,ifl)
424  zc = ppatch%fc(zcoord,ifl)
425 
426  fs = ppatch%gs(indgs*ifl)
427 
428  rl = cv(cv_mixt_dens,c1)
429  irl = 1.0_rfreal/rl
430 
431  ul = cv(cv_mixt_xmom,c1)*irl
432  vl = cv(cv_mixt_ymom,c1)*irl
433  wl = cv(cv_mixt_zmom,c1)*irl
434  tl = dv(dv_mixt_temp,c1)
435 
436  ygl = pcvspec(1,c1)*irl
437  yvl = pcvspec(2,c1)*irl
438 
439  dx = xc - pregion%grid%cofg(xcoord,c1)
440  dy = yc - pregion%grid%cofg(ycoord,c1)
441  dz = zc - pregion%grid%cofg(zcoord,c1)
442 
443  rl = rl + grad(xcoord,grc_mixt_dens,c1)*dx &
444  + grad(ycoord,grc_mixt_dens,c1)*dy &
445  + grad(zcoord,grc_mixt_dens,c1)*dz
446  ul = ul + grad(xcoord,grc_mixt_xvel,c1)*dx &
447  + grad(ycoord,grc_mixt_xvel,c1)*dy &
448  + grad(zcoord,grc_mixt_xvel,c1)*dz
449  vl = vl + grad(xcoord,grc_mixt_yvel,c1)*dx &
450  + grad(ycoord,grc_mixt_yvel,c1)*dy &
451  + grad(zcoord,grc_mixt_yvel,c1)*dz
452  wl = wl + grad(xcoord,grc_mixt_zvel,c1)*dx &
453  + grad(ycoord,grc_mixt_zvel,c1)*dy &
454  + grad(zcoord,grc_mixt_zvel,c1)*dz
455  tl = tl + grad(xcoord,grc_mixt_temp,c1)*dx &
456  + grad(ycoord,grc_mixt_temp,c1)*dy &
457  + grad(zcoord,grc_mixt_temp,c1)*dz
458 
459  ygl = ygl + pgradspec(xcoord,1,c1)*dx &
460  + pgradspec(ycoord,1,c1)*dy &
461  + pgradspec(zcoord,1,c1)*dz
462  yvl = yvl + pgradspec(xcoord,2,c1)*dx &
463  + pgradspec(ycoord,2,c1)*dy &
464  + pgradspec(zcoord,2,c1)*dz
465 
466  IF ( ygl > 1.0_rfreal ) THEN
467  ygl = 1.0_rfreal
468  yvl = 0.0_rfreal
469  ELSE IF ( ygl < 0.0_rfreal ) THEN
470  ygl = 0.0_rfreal
471  ELSE IF ( yvl > 1.0_rfreal ) THEN
472  yvl = 1.0_rfreal
473  ygl = 0.0_rfreal
474  ELSE IF ( yvl < 0.0_rfreal ) THEN
475  yvl = 0.0_rfreal
476  END IF ! Ygl
477 
478  yll = 1.0_rfreal - ygl - yvl
479 
480  IF ( yll > 1.0_rfreal) THEN
481  yll = 1.0_rfreal
482  ELSE IF ( yll < 0.0_rfreal) THEN
483  yll = 0.0_rfreal
484  END IF ! Yll
485 
486  vl2 = ul*ul + vl*vl + wl*wl
487 
488  cll2 = mixtliq_c2_bp(bp)
489  cvl2 = mixtperf_c2_grt(1.0_rfreal,rv,tl)
490  cgl2 = mixtperf_c2_grt(1.0_rfreal,rg,tl)
491 
492  rygl = rl*ygl
493  ryvl = rl*yvl
494  ryll = rl*yll
495 
496  pl = mixtgasliq_p(ryll,ryvl,rygl,cll2,cvl2,cgl2,rl,ro,po,to,bp,bt,tl)
497 
498  cvm = (ryll*cvl + ryvl*cvv + rygl*cvg)/rl
499 
500  rul = rl*ul
501  rvl = rl*vl
502  rwl = rl*wl
503  rel = rl*cvm*tl + 0.5*rl*vl2
504  rgpgl = rygl
505  rvpvl = ryvl
506  rlpll = ryll
507 
508 
509  IF ( bcopttype /= bcopt_supersonic ) THEN
510  pr = vals(bcdat_outflow_press,distrib*ifl)
511  ELSE
512  pr = pl
513  END IF ! bcOptType
514 
515  CALL bcondoutflowperf_gl(bcopttype,ro,po,to,bp,bt,cvl,cvv, &
516  cvg,rg,rv,pr,nx,ny,nz,rl,rul,rvl, &
517  rwl,rel,rgpgl,rvpvl,pl,rr,rur,rvr, &
518  rwr,rer,rgpgr,rvpvr)
519 
520  ql = (rul*nx + rvl*ny + rwl*nz)/rl - fs
521  qr = (rur*nx + rvr*ny + rwr*nz)/rr - fs
522 
523  flx(1) = 0.5_rfreal*(ql* rl + qr* rr )*nm
524  flx(2) = 0.5_rfreal*(ql* rul + pl*nx + qr* rur + pr*nx)*nm
525  flx(3) = 0.5_rfreal*(ql* rvl + pl*ny + qr* rvr + pr*ny)*nm
526  flx(4) = 0.5_rfreal*(ql* rwl + pl*nz + qr* rwr + pr*nz)*nm
527  flx(5) = 0.5_rfreal*(ql*(rel + pl) + pl*fs + qr*(rer + pr) + pr*fs)*nm
528 
529  mfmixt(indmf*ifl) = flx(1)
530  ppatch%cp(ifl) = icpref*(0.5_rfreal*(pl + pr) - pref)
531 
532  rhs(cv_mixt_dens,c1) = rhs(cv_mixt_dens,c1) + flx(1)
533  rhs(cv_mixt_xmom,c1) = rhs(cv_mixt_xmom,c1) + flx(2)
534  rhs(cv_mixt_ymom,c1) = rhs(cv_mixt_ymom,c1) + flx(3)
535  rhs(cv_mixt_zmom,c1) = rhs(cv_mixt_zmom,c1) + flx(4)
536  rhs(cv_mixt_ener,c1) = rhs(cv_mixt_ener,c1) + flx(5)
537 
538  sd(sd_xmom,c1*indsd) = sd(sd_xmom,c1*indsd) &
539  + 0.5_rfreal*(rul/rl + rur/rr)*flx(1)
540  sd(sd_ymom,c1*indsd) = sd(sd_ymom,c1*indsd) &
541  + 0.5_rfreal*(rvl/rl + rvr/rr)*flx(1)
542  sd(sd_zmom,c1*indsd) = sd(sd_zmom,c1*indsd) &
543  + 0.5_rfreal*(rwl/rl + rwr/rr)*flx(1)
544 
545  IF ( pregion%irkStep == 1 ) THEN
546  global%massIn = global%massIn - min(flx(1),0.0_rfreal)
547  global%massOut = global%massOut + max(flx(1),0.0_rfreal)
548  END IF ! pRegion
549 
550  IF ( (global%checkLevel == check_high) .AND. &
551  (global%verbLevel >= verbose_high) .AND. &
552  (global%myProcid == masterproc) .AND. &
553  (decideprintflag .EQV. .true.) ) THEN
554  IF ( flx(1) < 0.0_rfreal ) THEN
555  nlocs = nlocs + 1
556 
557  IF ( nlocs == 1 ) THEN
558  global%warnCounter = global%warnCounter + 1
559 
560  WRITE(stdout,'(A,1X,A,1X,I9)') solver_name, &
561  '*** WARNING *** Inflow detected at outflow boundary!'
562  WRITE(stdout,'(A,3X,A,1X,I5.5)') solver_name,'Global region:', &
563  pregion%iRegionGlobal
564  IF ( global%flowType == flow_unsteady ) THEN
565  WRITE(stdout,'(A,3X,A,1X,1PE11.5)') solver_name, &
566  'Current time:', &
567  global%currentTime
568  ELSE
569  WRITE(stdout,'(A,3X,A,1X,I6.6)') solver_name, &
570  'Current iteration number:', &
571  global%currentIter
572  END IF ! global%flowType
573 
574  WRITE(stdout,'(A,3X,A,1X,I1)') solver_name, &
575  'Runge-Kutta stage:', &
576  pregion%irkStep
577  END IF ! nLocs
578 
579  IF ( nlocs <= max_inoutflow_locs ) THEN
580  loc(nlocs,min_val:max_val) = c1
581  END IF ! nLocs
582  END IF ! flx(1)
583  END IF ! global%checkLevel
584  END DO ! ifl
585 
586 ! ------------------------------------------------------------------------------
587 ! Write info on inflow at outflow boundary
588 ! ------------------------------------------------------------------------------
589 
590  IF ( (global%checkLevel == check_high) .AND. &
591  (global%verbLevel >= verbose_high) .AND. &
592  (global%myProcid == masterproc) .AND. &
593  (decideprintflag .EQV. .true.) .AND. &
594  (nlocs > 0) ) THEN
595  IF ( nlocs > max_inoutflow_locs ) THEN
596  WRITE(stdout,'(A,3X,A,1X,I3,1X,A,1X,I9,1X,A)') solver_name, &
597  'Only wrote the first',max_inoutflow_locs,'of',nlocs, &
598  'outflow faces with inflow.'
599  CALL rflu_printlocinfo(pregion,loc,max_inoutflow_locs, &
600  locinfo_mode_silent,output_mode_anybody)
601  ELSE
602  CALL rflu_printlocinfo(pregion,loc(1:nlocs,min_val:max_val),nlocs, &
603  locinfo_mode_silent,output_mode_anybody)
604  END IF ! nLocs
605  END IF ! global%verbLevel
606 
607 ! ==============================================================================
608 ! Slip wall (weak imposition)
609 ! ==============================================================================
610 
611  CASE ( bc_slipwall )
612  DO ifl = 1,ppatch%nBFaces
613  c1 = ppatch%bf2c(ifl)
614 
615  nx = ppatch%fn(xcoord,ifl)
616  ny = ppatch%fn(ycoord,ifl)
617  nz = ppatch%fn(zcoord,ifl)
618  nm = ppatch%fn(xyzmag,ifl)
619 
620  xc = ppatch%fc(xcoord,ifl)
621  yc = ppatch%fc(ycoord,ifl)
622  zc = ppatch%fc(zcoord,ifl)
623 
624  fs = ppatch%gs(indgs*ifl)
625  fsu = rflu_descalegridspeed(pregion,fs)
626 
627  rl = cv(cv_mixt_dens,c1)
628  irl = 1.0_rfreal/rl
629 
630  ul = cv(cv_mixt_xmom,c1)*irl
631  vl = cv(cv_mixt_ymom,c1)*irl
632  wl = cv(cv_mixt_zmom,c1)*irl
633  tl = dv(dv_mixt_temp,c1)
634 
635  ygl = pcvspec(1,c1)*irl
636  yvl = pcvspec(2,c1)*irl
637 
638  dx = xc - pregion%grid%cofg(xcoord,c1)
639  dy = yc - pregion%grid%cofg(ycoord,c1)
640  dz = zc - pregion%grid%cofg(zcoord,c1)
641 
642  rl = rl + grad(xcoord,grc_mixt_dens,c1)*dx &
643  + grad(ycoord,grc_mixt_dens,c1)*dy &
644  + grad(zcoord,grc_mixt_dens,c1)*dz
645  ul = ul + grad(xcoord,grc_mixt_xvel,c1)*dx &
646  + grad(ycoord,grc_mixt_xvel,c1)*dy &
647  + grad(zcoord,grc_mixt_xvel,c1)*dz
648  vl = vl + grad(xcoord,grc_mixt_yvel,c1)*dx &
649  + grad(ycoord,grc_mixt_yvel,c1)*dy &
650  + grad(zcoord,grc_mixt_yvel,c1)*dz
651  wl = wl + grad(xcoord,grc_mixt_zvel,c1)*dx &
652  + grad(ycoord,grc_mixt_zvel,c1)*dy &
653  + grad(zcoord,grc_mixt_zvel,c1)*dz
654  tl = tl + grad(xcoord,grc_mixt_temp,c1)*dx &
655  + grad(ycoord,grc_mixt_temp,c1)*dy &
656  + grad(zcoord,grc_mixt_temp,c1)*dz
657 
658  ygl = ygl + pgradspec(xcoord,1,c1)*dx &
659  + pgradspec(ycoord,1,c1)*dy &
660  + pgradspec(zcoord,1,c1)*dz
661  yvl = yvl + pgradspec(xcoord,2,c1)*dx &
662  + pgradspec(ycoord,2,c1)*dy &
663  + pgradspec(zcoord,2,c1)*dz
664 
665  IF ( ygl > 1.0_rfreal ) THEN
666  ygl = 1.0_rfreal
667  yvl = 0.0_rfreal
668  ELSE IF ( ygl < 0.0_rfreal ) THEN
669  ygl = 0.0_rfreal
670  ELSE IF ( yvl > 1.0_rfreal ) THEN
671  yvl = 1.0_rfreal
672  ygl = 0.0_rfreal
673  ELSE IF ( yvl < 0.0_rfreal ) THEN
674  yvl = 0.0_rfreal
675  END IF ! Ygl
676 
677  yll = 1.0_rfreal - ygl - yvl
678 
679  IF ( yll > 1.0_rfreal) THEN
680  yll = 1.0_rfreal
681  ELSE IF ( yll < 0.0_rfreal) THEN
682  yll = 0.0_rfreal
683  END IF ! Yll
684 
685  cll2 = mixtliq_c2_bp(bp)
686  cvl2 = mixtperf_c2_grt(1.0_rfreal,rv,tl)
687  cgl2 = mixtperf_c2_grt(1.0_rfreal,rg,tl)
688 
689  rygl = rl*ygl
690  ryvl = rl*yvl
691  ryll = rl*yll
692 
693  pl = mixtgasliq_p(ryll,ryvl,rygl,cll2,cvl2,cgl2,rl,ro,po,to,bp,bt,tl)
694 
695  cvm = (ryll*cvl + ryvl*cvv + rygl*cvg)/rl
696 
697  rul = rl*ul
698  rvl = rl*vl
699  rwl = rl*wl
700 
701  cp = gv(gv_mixt_cp ,indcp *c1)
702  mm = gv(gv_mixt_mol,indmol*c1)
703 
704  CALL rflu_setrindstateslipwallperf(cp,mm,nx,ny,nz,rl,rul, &
705  rvl,rwl,fsu,pl)
706 
707  flx(2) = pl*nx*nm
708  flx(3) = pl*ny*nm
709  flx(4) = pl*nz*nm
710  flx(5) = pl*fs*nm
711 
712  mfmixt(indmf*ifl) = 0.0_rfreal
713  ppatch%cp(ifl) = icpref*(pl - pref)
714 
715  rhs(cv_mixt_xmom,c1) = rhs(cv_mixt_xmom,c1) + flx(2)
716  rhs(cv_mixt_ymom,c1) = rhs(cv_mixt_ymom,c1) + flx(3)
717  rhs(cv_mixt_zmom,c1) = rhs(cv_mixt_zmom,c1) + flx(4)
718  rhs(cv_mixt_ener,c1) = rhs(cv_mixt_ener,c1) + flx(5)
719  END DO ! ifl
720 
721 ! ==============================================================================
722 ! No-slip wall
723 ! ==============================================================================
724 
725  CASE ( bc_noslipwall )
726  DO ifl = 1,ppatch%nBFaces
727  c1 = ppatch%bf2c(ifl)
728 
729  nx = ppatch%fn(xcoord,ifl)
730  ny = ppatch%fn(ycoord,ifl)
731  nz = ppatch%fn(zcoord,ifl)
732  nm = ppatch%fn(xyzmag,ifl)
733 
734  xc = ppatch%fc(xcoord,ifl)
735  yc = ppatch%fc(ycoord,ifl)
736  zc = ppatch%fc(zcoord,ifl)
737 
738  fs = ppatch%gs(indgs*ifl)
739 
740  rl = cv(cv_mixt_dens,c1)
741  tl = dv(dv_mixt_temp,c1)
742  ygl = pcvspec(1,c1)*irl
743  yvl = pcvspec(2,c1)*irl
744 
745  dx = xc - pregion%grid%cofg(xcoord,c1)
746  dy = yc - pregion%grid%cofg(ycoord,c1)
747  dz = zc - pregion%grid%cofg(zcoord,c1)
748 
749  rl = rl + grad(xcoord,grc_mixt_dens,c1)*dx &
750  + grad(ycoord,grc_mixt_dens,c1)*dy &
751  + grad(zcoord,grc_mixt_dens,c1)*dz
752  tl = tl + grad(xcoord,grc_mixt_temp,c1)*dx &
753  + grad(ycoord,grc_mixt_temp,c1)*dy &
754  + grad(zcoord,grc_mixt_temp,c1)*dz
755 
756  ygl = ygl + pgradspec(xcoord,1,c1)*dx &
757  + pgradspec(ycoord,1,c1)*dy &
758  + pgradspec(zcoord,1,c1)*dz
759  yvl = yvl + pgradspec(xcoord,2,c1)*dx &
760  + pgradspec(ycoord,2,c1)*dy &
761  + pgradspec(zcoord,2,c1)*dz
762 
763  IF ( ygl > 1.0_rfreal ) THEN
764  ygl = 1.0_rfreal
765  yvl = 0.0_rfreal
766  ELSE IF ( ygl < 0.0_rfreal ) THEN
767  ygl = 0.0_rfreal
768  ELSE IF ( yvl > 1.0_rfreal ) THEN
769  yvl = 1.0_rfreal
770  ygl = 0.0_rfreal
771  ELSE IF ( yvl < 0.0_rfreal ) THEN
772  yvl = 0.0_rfreal
773  END IF ! Ygl
774 
775  yll = 1.0_rfreal - ygl - yvl
776 
777  IF ( yll > 1.0_rfreal) THEN
778  yll = 1.0_rfreal
779  ELSE IF ( yll < 0.0_rfreal) THEN
780  yll = 0.0_rfreal
781  END IF ! Yll
782 
783  cll2 = mixtliq_c2_bp(bp)
784  cvl2 = mixtperf_c2_grt(1.0_rfreal,rv,tl)
785  cgl2 = mixtperf_c2_grt(1.0_rfreal,rg,tl)
786 
787  rygl = rl*ygl
788  ryvl = rl*yvl
789  ryll = rl*yll
790 
791  pl = mixtgasliq_p(ryll,ryvl,rygl,cll2,cvl2,cgl2,rl,ro,po,to,bp,bt,tl)
792 
793  flx(2) = pl*nx*nm
794  flx(3) = pl*ny*nm
795  flx(4) = pl*nz*nm
796  flx(5) = pl*fs*nm
797 
798  mfmixt(indmf*ifl) = 0.0_rfreal
799  ppatch%cp(ifl) = icpref*(pl - pref)
800 
801  rhs(cv_mixt_xmom,c1) = rhs(cv_mixt_xmom,c1) + flx(2)
802  rhs(cv_mixt_ymom,c1) = rhs(cv_mixt_ymom,c1) + flx(3)
803  rhs(cv_mixt_zmom,c1) = rhs(cv_mixt_zmom,c1) + flx(4)
804  rhs(cv_mixt_ener,c1) = rhs(cv_mixt_ener,c1) + flx(5)
805  END DO ! ifl
806 
807 ! ==============================================================================
808 ! Boundaries for which fluxes must not or need not be computed
809 ! ==============================================================================
810 
811  CASE ( bc_periodic, &
812  bc_symmetry, &
813  bc_virtual )
814 
815 ! ==============================================================================
816 ! Default
817 ! ==============================================================================
818 
819  CASE default
820  CALL errorstop(global,err_reached_default,__line__)
821  END SELECT ! pPatch%bcType
822 
823 ! ******************************************************************************
824 ! End
825 ! ******************************************************************************
826 
827  CALL deregisterfunction(global)
828 
829 END SUBROUTINE rflu_centralsecondpatch_gl
830 
831 ! ******************************************************************************
832 !
833 ! RCS Revision history:
834 !
835 ! $Log: RFLU_CentralSecondPatch_GL.F90,v $
836 ! Revision 1.5 2008/12/06 08:44:29 mtcampbe
837 ! Updated license.
838 !
839 ! Revision 1.4 2008/11/19 22:17:42 mtcampbe
840 ! Added Illinois Open Source License/Copyright
841 !
842 ! Revision 1.3 2006/08/19 15:39:56 mparmar
843 ! Renamed patch variables
844 !
845 ! Revision 1.2 2006/04/13 20:02:27 haselbac
846 ! Added periodic and virtual boundary CASE statement
847 !
848 ! Revision 1.1 2006/03/26 20:21:03 haselbac
849 ! Initial revision
850 !
851 ! ******************************************************************************
852 
853 
854 
855 
856 
857 
858 
subroutine, public rflu_setrindstateinjectperf(cpGas, mmGas, nx, ny, nz, mInj, tInj, pl, fs, rl, ul, vl, wl, Hl)
real(rfreal) function mixtperf_p_deogvm2(D, Eo, G, Vm2)
Definition: MixtPerf_P.F90:39
real(rfreal) function mixtperf_r_m(M)
Definition: MixtPerf_R.F90:54
NT dx
NT rhs
subroutine, public rflu_setrindstateslipwallperf(cpGas, mmGas, nx, ny, nz, rl, rul, rvl, rwl, fs, pl)
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
real(rfreal) function mixtperf_d_prt(P, R, T)
Definition: MixtPerf_D.F90:71
double sqrt(double d)
Definition: double.h:73
real(rfreal) function mixtperf_r_cpg(Cp, G)
Definition: MixtPerf_R.F90:39
real(rfreal) function mixtgasliq_p(DYl, DYv, DYg, Cl2, Cv2, Cg2, D, Dz, Po, To, Bp, Bt, T)
subroutine rflu_centralsecondpatch_gl(pRegion, pPatch)
real(rfreal) function, public rflu_descalegridspeed(pRegion, fs)
real(rfreal) function mixtperf_c2_grt(G, R, T)
Definition: MixtPerf_C.F90:101
real(rfreal) function mixtliq_c2_bp(Bp)
Definition: MixtLiq_C.F90:54
subroutine, public rflu_setrindstatefarfieldperf(global, cpGas, mmGas, nx, ny, nz, machInf, pInf, tInf, alphaInf, betaInf, corrFlag, liftCoef, xc, yc, zc, rl, rul, rvl, rwl, rel, rr, rur, rvr, rwr, rer, pr)
RT dz() const
Definition: Direction_3.h:133
Vector_n min(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:346
real(rfreal) function mixtliq_d_dobpppobttto(Dz, Bp, Bt, P, Po, T, To)
Definition: MixtLiq_D.F90:40
subroutine rflu_printlocinfo(pRegion, locUnsorted, nLocUnsorted, locInfoMode, outputMode)
subroutine bcondinflowperf_gl(bcOptType, ro, po, to, Bp, Bt, cvl, cvv, cvg, Rg, Rv, ur, vr, wr, vfgr, vfvr, vflr, temp, press, nx, ny, nz, rl, rul, rvl, rwl, rel, rgpgl, rvpvl, pl, rr, rur, rvr, rwr, rer, rgpgr, rvpvr, pr)
real(rfreal) function mixtgasliq_c(Cvm, D, P, Dl, Dv, Dg, VFl, VFv, VFg, Cl2, Cv2, Cg2, Bl2, Bv2, Bg2)
NT dy
LOGICAL function rflu_decideprint(global)
real(rfreal) function mixtperf_eo_dgpuvw(D, G, P, U, V, W)
Definition: MixtPerf_E.F90:40
subroutine errorstop(global, errorCode, errorLine, addMessage)
Definition: ModError.F90:483
subroutine bcondoutflowperf_gl(bcOpt, ro, Po, To, betaP, betaT, cvl, cvv, cvg, Rg, Rv, pout, sxn, syn, szn, rho, rhou, rhov, rhow, rhoe, rhogpg, rhovpv, pin, rhob, rhoub, rhovb, rhowb, rhoeb, rhogpgb, rhovpvb)
subroutine deregisterfunction(global)
Definition: ModError.F90:469
real(rfreal) function mixtperf_g_cpr(Cp, R)
Definition: MixtPerf_G.F90:39
real(rfreal) function mixtperf_cv_cpr(Cp, R)
Definition: MixtPerf_Cv.F90:39
real(rfreal) function mixtperf_p_drt(D, R, T)
Definition: MixtPerf_P.F90:54