Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLU_ModHLLCFlux.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: Collection of HLLC flux routines.
26 !
27 ! Description: None.
28 !
29 ! Notes: None.
30 !
31 ! ******************************************************************************
32 !
33 ! $Id: RFLU_ModHLLCFlux.F90,v 1.10 2008/12/06 08:44:22 mtcampbe Exp $
34 !
35 ! Copyright: (c) 2005-2006 by the University of Illinois
36 !
37 ! ******************************************************************************
38 
39 MODULE rflu_modhllcflux
40 
41  USE moddatatypes
42  USE modglobal, ONLY: t_global
43  USE modparameters
44  USE moderror
45  USE moddatastruct, ONLY: t_region
46  USE modbndpatch, ONLY: t_patch
47  USE modgrid, ONLY: t_grid
48  USE modmpi
49 
50  IMPLICIT NONE
51 
52  PRIVATE
53  PUBLIC :: rflu_hllc_computeflux
54 
55 ! ******************************************************************************
56 ! Declarations and definitions
57 ! ******************************************************************************
58 
59  CHARACTER(CHRLEN) :: &
60  RCSIdentString = '$RCSfile: RFLU_ModHLLCFlux.F90,v $ $Revision: 1.10 $'
61 
62 ! ******************************************************************************
63 ! Routines
64 ! ******************************************************************************
65 
66  CONTAINS
67 
68 
69 
70 
71 
72 
73 
74 ! ******************************************************************************
75 !
76 ! Purpose: Wrapper function for HLLC flux functions.
77 !
78 ! Description: None.
79 !
80 ! Input:
81 ! pRegion Pointer to region data.
82 !
83 ! Output: None.
84 !
85 ! Notes: None.
86 !
87 ! ******************************************************************************
88 
89 SUBROUTINE rflu_hllc_computeflux(pRegion)
90 
91  IMPLICIT NONE
92 
93 ! ******************************************************************************
94 ! Definitions and declarations
95 ! ******************************************************************************
96 
97 ! ==============================================================================
98 ! Arguments
99 ! ==============================================================================
100 
101  TYPE(t_region), POINTER :: pregion
102 
103 ! ==============================================================================
104 ! Locals
105 ! ==============================================================================
106 
107  TYPE(t_global), POINTER :: global
108 
109 ! ******************************************************************************
110 ! Start
111 ! ******************************************************************************
112 
113  global => pregion%global
114 
115  CALL registerfunction(global,'RFLU_HLLC_ComputeFlux',&
116  'RFLU_ModHLLCFlux.F90')
117 
118 ! ******************************************************************************
119 ! Call flux functions
120 ! ******************************************************************************
121 
122  SELECT CASE ( pregion%mixtInput%gasModel )
123 
124 ! ==============================================================================
125 ! Thermally and calorically perfect gas
126 ! ==============================================================================
127 
128  CASE ( gas_model_tcperf )
129  SELECT CASE ( pregion%mixtInput%spaceOrder )
130  CASE ( discr_order_1 )
131  CALL rflu_hllc_computeflux1_tcp(pregion)
132  CASE ( discr_order_2 )
133  CALL rflu_hllc_computeflux2_tcp(pregion)
134  CASE default
135  CALL errorstop(global,err_reached_default,__line__)
136  END SELECT ! pRegion%mixtInput%spaceOrder
137 
138 ! ==============================================================================
139 ! Mixture of thermally and calorically perfect gases
140 ! ==============================================================================
141 
142  CASE ( gas_model_mixt_tcperf )
143  SELECT CASE ( pregion%mixtInput%spaceOrder )
144  CASE ( discr_order_1 )
145  CALL rflu_hllc_computeflux1_mtcp(pregion)
146  CASE ( discr_order_2 )
147  CALL rflu_hllc_computeflux2_mtcp(pregion)
148  CASE default
149  CALL errorstop(global,err_reached_default,__line__)
150  END SELECT ! pRegion%mixtInput%spaceOrder
151 
152 ! ==============================================================================
153 ! Mixture of gas, vapor, and liquid
154 ! ==============================================================================
155 
156  CASE ( gas_model_mixt_gasliq )
157  SELECT CASE ( pregion%mixtInput%spaceOrder )
158  CASE ( discr_order_1 )
159  CALL rflu_hllc_computeflux1_gl(pregion)
160  CASE ( discr_order_2 )
161  CALL rflu_hllc_computeflux2_gl(pregion)
162  CASE default
163  CALL errorstop(global,err_reached_default,__line__)
164  END SELECT ! pRegion%mixtInput%spaceOrder
165 
166 ! ==============================================================================
167 ! Default
168 ! ==============================================================================
169 
170  CASE default
171  CALL errorstop(global,err_reached_default,__line__)
172  END SELECT ! pRegion%mixtInput%gasModel
173 
174 ! ******************************************************************************
175 ! End
176 ! ******************************************************************************
177 
178  CALL deregisterfunction(global)
179 
180 END SUBROUTINE rflu_hllc_computeflux
181 
182 
183 
184 
185 
186 ! ******************************************************************************
187 !
188 ! Purpose: Compute central convective fluxes for multiphase mixture using
189 ! first-order accurate version of HLLC scheme.
190 !
191 ! Description: None.
192 !
193 ! Input:
194 ! pRegion Data of current region.
195 !
196 ! Output: None.
197 !
198 ! Notes: None.
199 !
200 ! ******************************************************************************
201 
202 SUBROUTINE rflu_hllc_computeflux1_gl(pRegion)
203 
204  USE modinterfaces, ONLY: mixtgasliq_c, &
205  mixtliq_d_dobpppobttto, &
206  mixtliq_c2_bp, &
207  mixtperf_c2_grt, &
208  mixtperf_cv_cpr, &
209  mixtperf_d_prt, &
210  mixtperf_r_cpg, &
211  mixtperf_r_m, &
212  rflu_centralfirstpatch_gl
213 
214  IMPLICIT NONE
215 
216 ! ******************************************************************************
217 ! Definitions and declarations
218 ! ******************************************************************************
219 
220 ! ==============================================================================
221 ! Arguments
222 ! ==============================================================================
223 
224  TYPE(t_region), POINTER :: pregion
225 
226 ! ==============================================================================
227 ! Locals
228 ! ==============================================================================
229 
230  INTEGER :: c1,c2,ifg,indgs,indmf,indsd,ipatch
231  INTEGER, DIMENSION(:,:), POINTER :: f2c
232  REAL(RFREAL) :: bgh2,blh2,bp,bt,bvh2,cml,cmr,cms,cvg,cvl,cvv,cgh2,clh2, &
233  cvh2,cvm,el,er,fs,irl,irr,nm,nx,ny,nz,ph,pl,po,pr,ps,ql, &
234  qr,qs,res,rgh,rgl,rgr,rh,rl,rlh,ro,rr,rs,rus,rvh,rvl,rvr, &
235  rvs,rws,rygh,rygl,rygr,ryvh,ryvl,ryvr,rg,rv,sl,sm,sr, &
236  term,th,tl,to,tr,ul,ur,us,vfgh,vfgl,vfgr,vflh,vfvh,vfvl, &
237  vfvr,vl,vr,vs,vms,wl,wr,ws,wt1,wt2
238  REAL(RFREAL) :: flx(5)
239  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
240  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcvmixt,pcvspec,pdvmixt,prhs,psd
241  TYPE(t_global), POINTER :: global
242  TYPE(t_grid), POINTER :: pgrid
243  TYPE(t_patch), POINTER :: ppatch
244 
245 ! ******************************************************************************
246 ! Start
247 ! ******************************************************************************
248 
249  global => pregion%global
250 
251  CALL registerfunction(global,'RFLU_HLLC_ComputeFlux1_GL',&
252  'RFLU_ModHLLCFlux.F90')
253 
254 ! ******************************************************************************
255 ! Set dimensions and pointers
256 ! ******************************************************************************
257 
258  pgrid => pregion%grid
259 
260  indgs = pgrid%indGs
261 
262  indmf = pregion%mixtInput%indMfMixt
263  indsd = pregion%mixtInput%indSd
264 
265  pcvmixt => pregion%mixt%cv
266  pdvmixt => pregion%mixt%dv
267  pcvspec => pregion%spec%cv
268 
269  pmf => pregion%mixt%mfMixt
270  prhs => pregion%mixt%rhs
271  psd => pregion%mixt%sd
272 
273  IF ( pregion%spec%cvState /= cv_mixt_state_cons ) THEN
274  CALL errorstop(global,err_cv_state_invalid,__line__)
275  END IF ! pRegion%spec%cvState
276 
277  IF ( pregion%mixtInput%gasModel /= gas_model_mixt_gasliq ) THEN
278  CALL errorstop(global,err_gasmodel_invalid,__line__)
279  END IF ! pRegion%mixtInput%gasModel
280 
281 ! ******************************************************************************
282 ! Define constants
283 ! ******************************************************************************
284 
285  ro = global%refDensityLiq
286  po = global%refPressLiq
287  to = global%refTempLiq
288  bp = global%refBetaPLiq
289  bt = global%refBetaTLiq
290  cvl = global%refCvLiq
291 
292  rg = mixtperf_r_m(pregion%specInput%specType(1)%pMaterial%molw)
293  cvg = mixtperf_cv_cpr(pregion%specInput%specType(1)%pMaterial%spht,rg)
294 
295  rv = mixtperf_r_m(pregion%specInput%specType(2)%pMaterial%molw)
296  cvv = mixtperf_cv_cpr(pregion%specInput%specType(2)%pMaterial%spht,rv)
297 
298 ! ******************************************************************************
299 ! Compute fluxes through interior faces
300 ! ******************************************************************************
301 
302  DO ifg = 1,pgrid%nFaces
303  c1 = pgrid%f2c(1,ifg)
304  c2 = pgrid%f2c(2,ifg)
305 
306 ! ==============================================================================
307 ! Get face geometry and grid speed
308 ! ==============================================================================
309 
310  nx = pgrid%fn(xcoord,ifg)
311  ny = pgrid%fn(ycoord,ifg)
312  nz = pgrid%fn(zcoord,ifg)
313  nm = pgrid%fn(xyzmag,ifg)
314 
315  fs = pgrid%gs(indgs*ifg)
316 
317 ! =============================================================================
318 ! Compute left and right states
319 ! =============================================================================
320 
321 ! -----------------------------------------------------------------------------
322 ! Left state
323 ! -----------------------------------------------------------------------------
324 
325  rl = pcvmixt(cv_mixt_dens,c1)
326  irl = 1.0_rfreal/rl
327 
328  ul = pcvmixt(cv_mixt_xmom,c1)*irl
329  vl = pcvmixt(cv_mixt_ymom,c1)*irl
330  wl = pcvmixt(cv_mixt_zmom,c1)*irl
331  el = pcvmixt(cv_mixt_ener,c1)*irl
332  pl = pdvmixt(dv_mixt_pres,c1)
333  tl = pdvmixt(dv_mixt_temp,c1)
334  cml = pdvmixt(dv_mixt_soun,c1)
335 
336  rygl = pcvspec(1,c1)
337  ryvl = pcvspec(2,c1)
338 
339  rvl = mixtperf_d_prt(pl,rv,tl)
340  rgl = mixtperf_d_prt(pl,rg,tl)
341 
342  vfvl = ryvl/rvl
343  vfgl = rygl/rgl
344 
345  ql = ul*nx + vl*ny + wl*nz - fs
346 
347 ! -----------------------------------------------------------------------------
348 ! Right state
349 ! -----------------------------------------------------------------------------
350 
351  rr = pcvmixt(cv_mixt_dens,c2)
352  irr = 1.0_rfreal/rr
353 
354  ur = pcvmixt(cv_mixt_xmom,c2)*irr
355  vr = pcvmixt(cv_mixt_ymom,c2)*irr
356  wr = pcvmixt(cv_mixt_zmom,c2)*irr
357  er = pcvmixt(cv_mixt_ener,c2)*irr
358  pr = pdvmixt(dv_mixt_pres,c2)
359  tr = pdvmixt(dv_mixt_temp,c2)
360  cmr = pdvmixt(dv_mixt_soun,c2)
361 
362  rygr = pcvspec(1,c2)
363  ryvr = pcvspec(2,c2)
364 
365  rvr = mixtperf_d_prt(pr,rv,tr)
366  rgr = mixtperf_d_prt(pr,rg,tr)
367 
368  vfvr = ryvr/rvr
369  vfgr = rygr/rgr
370 
371  qr = ur*nx + vr*ny + wr*nz - fs
372 
373 ! ==============================================================================
374 ! Compute weights
375 ! ==============================================================================
376 
377  wt1 = sqrt(rr/rl)
378  wt2 = 1.0_rfreal/(1.0_rfreal + wt1)
379 
380 ! ==============================================================================
381 ! Compute averages
382 ! ==============================================================================
383 
384  us = wt2*(ul + wt1*ur)
385  vs = wt2*(vl + wt1*vr)
386  ws = wt2*(wl + wt1*wr)
387  qs = wt2*(ql + wt1*qr)
388 
389  vms = us*us + vs*vs + ws*ws
390 
391  rh = wt1*rl
392  ph = wt2*(pl + wt1*pr)
393  th = wt2*(tl + wt1*tr)
394  rygh = wt2*(rygl + wt1*rygr)
395  ryvh = wt2*(ryvl + wt1*ryvr)
396 
397  rlh = mixtliq_d_dobpppobttto(ro,bp,bt,ph,po,th,to)
398  rvh = mixtperf_d_prt(ph,rv,th)
399  rgh = mixtperf_d_prt(ph,rg,th)
400 
401  vfgh = rygh/rgh
402  vfvh = ryvh/rvh
403 
404  IF ( vfgh > 1.0_rfreal) THEN
405  vfgh = 1.0_rfreal
406  vfvh = 0.0_rfreal
407  ELSE IF ( vfgh < 0.0_rfreal) THEN
408  vfgh = 0.0_rfreal
409  ELSE IF ( vfvh > 1.0_rfreal) THEN
410  vfvh = 1.0_rfreal
411  vfgh = 0.0_rfreal
412  ELSE IF ( vfvh < 0.0_rfreal) THEN
413  vfvh = 0.0_rfreal
414  END IF ! vfgh
415 
416  vflh = 1.0_rfreal - vfvh - vfgh
417 
418  IF ( vflh > 1.0_rfreal ) THEN
419  vflh = 1.0_rfreal
420  ELSE IF ( vflh < 0.0_rfreal) THEN
421  vflh = 0.0_rfreal
422  END IF ! vflh
423 
424  clh2 = mixtliq_c2_bp(bp)
425  cvh2 = mixtperf_c2_grt(1.0_rfreal,rv,th)
426  cgh2 = mixtperf_c2_grt(1.0_rfreal,rg,th)
427 
428  blh2 = -bt/bp
429  bvh2 = rvh*rv
430  bgh2 = rgh*rg
431 
432  cvm = (rlh*vflh*cvl + rvh*vfvh*cvv + rgh*vfgh*cvg)/rh
433 
434  cms = mixtgasliq_c(cvm,rh,ph,rlh,rvh,rgh,vflh,vfvh,vfgh,clh2,cvh2,cgh2, &
435  blh2,bvh2,bgh2)
436 
437  sl = min(ql-cml,qs-cms)
438  sr = max(qr+cmr,qs+cms)
439 
440 ! ==============================================================================
441 ! Compute fluxes
442 ! ==============================================================================
443 
444  IF ( sl > 0.0_rfreal ) THEN
445  flx(1) = ql* rl *nm
446  flx(2) = (ql* rl*ul + pl*nx )*nm
447  flx(3) = (ql* rl*vl + pl*ny )*nm
448  flx(4) = (ql* rl*wl + pl*nz )*nm
449  flx(5) = (ql*(rl*el + pl) + pl*fs)*nm
450  ELSE IF ( sr < 0.0_rfreal ) THEN
451  flx(1) = qr* rr *nm
452  flx(2) = (qr* rr*ur + pr*nx )*nm
453  flx(3) = (qr* rr*vr + pr*ny )*nm
454  flx(4) = (qr* rr*wr + pr*nz )*nm
455  flx(5) = (qr*(rr*er + pr) + pr*fs)*nm
456  ELSE
457  sm = (rr*qr*(sr-qr) - rl*ql*(sl-ql) - pr + pl)/(rr*(sr-qr) - rl*(sl-ql))
458  ps = pl + rl*(ql-sl)*(ql-sm)
459 
460  IF ( sm > 0.0_rfreal ) THEN
461  term = 1.0_rfreal/(sl-sm)
462 
463  rs = term* (sl-ql)*rl
464  rus = term*((sl-ql)*rl*ul + (ps - pl )*nx)
465  rvs = term*((sl-ql)*rl*vl + (ps - pl )*ny)
466  rws = term*((sl-ql)*rl*wl + (ps - pl )*nz)
467  res = term*((sl-ql)*rl*el + (ps*sm - pl*ql) )
468  ELSE
469  term = 1.0_rfreal/(sr-sm)
470 
471  rs = term* (sr-qr)*rr
472  rus = term*((sr-qr)*rr*ur + (ps - pr )*nx)
473  rvs = term*((sr-qr)*rr*vr + (ps - pr )*ny)
474  rws = term*((sr-qr)*rr*wr + (ps - pr )*nz)
475  res = term*((sr-qr)*rr*er + (ps*sm - pr*qr) )
476  END IF ! sm
477 
478  flx(1) = sm* rs *nm
479  flx(2) = (sm* rus + ps*nx )*nm
480  flx(3) = (sm* rvs + ps*ny )*nm
481  flx(4) = (sm* rws + ps*nz )*nm
482  flx(5) = (sm*(res + ps) + ps*fs)*nm
483  END IF ! sl
484 
485 ! ==============================================================================
486 ! Store mass flux
487 ! ==============================================================================
488 
489  pmf(indmf*ifg) = flx(1)
490 
491 ! ==============================================================================
492 ! Accumulate into residual
493 ! ==============================================================================
494 
495  prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + flx(1)
496  prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + flx(2)
497  prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + flx(3)
498  prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + flx(4)
499  prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + flx(5)
500 
501  prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - flx(1)
502  prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - flx(2)
503  prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - flx(3)
504  prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - flx(4)
505  prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - flx(5)
506 
507 ! ==============================================================================
508 ! Accumulate into substantial derivative
509 ! ==============================================================================
510 
511  psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + us*pmf(indmf*ifg)
512  psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vs*pmf(indmf*ifg)
513  psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + ws*pmf(indmf*ifg)
514 
515  psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - us*pmf(indmf*ifg)
516  psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vs*pmf(indmf*ifg)
517  psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - ws*pmf(indmf*ifg)
518  END DO ! ifg
519 
520 ! ******************************************************************************
521 ! Compute fluxes through boundary faces
522 ! ******************************************************************************
523 
524  DO ipatch = 1,pregion%grid%nPatches
525  ppatch => pregion%patches(ipatch)
526 
527  CALL rflu_centralfirstpatch_gl(pregion,ppatch)
528  END DO ! iPatch
529 
530 ! ******************************************************************************
531 ! End
532 ! ******************************************************************************
533 
534  CALL deregisterfunction(global)
535 
536 END SUBROUTINE rflu_hllc_computeflux1_gl
537 
538 
539 
540 
541 
542 
543 
544 
545 ! ******************************************************************************
546 !
547 ! Purpose: Compute convective fluxes using first-order accurate version of
548 ! HLLC scheme for mixture of thermally and calorically perfect gases.
549 !
550 ! Description: None.
551 !
552 ! Input:
553 ! pRegion Pointer to region
554 !
555 ! Output: None.
556 !
557 ! Notes: None.
558 !
559 ! ******************************************************************************
560 
561 SUBROUTINE rflu_hllc_computeflux1_mtcp(pRegion)
562 
563 #ifdef SPEC
564  USE modspecies, ONLY: t_spec_type
565 #endif
566 
567  USE modinterfaces, ONLY: mixtperf_c_dgp, &
568  mixtperf_c_ghovm2, &
569  mixtperf_eo_dgpuvw, &
570  mixtperf_g_cpr, &
571  mixtperf_r_m, &
572  rflu_centralfirstpatch
573 
574  IMPLICIT NONE
575 
576 ! ******************************************************************************
577 ! Definitions and declarations
578 ! ******************************************************************************
579 
580 ! ==============================================================================
581 ! Arguments
582 ! ==============================================================================
583 
584  TYPE(t_region), POINTER :: pregion
585 
586 ! ==============================================================================
587 ! Locals
588 ! ==============================================================================
589 
590  INTEGER :: c1,c2,ifg,indcp,indgs,indmf,indmol,indsd,ipatch
591  REAL(RFREAL) :: al,ar,as,cpl,cpr,el,er,fs,gcl,gcr,gl,gr,gs,hl,hr,hs,irl, &
592  irr,mml,mmr,nm,nx,ny,nz,ql,qr,qs,pl,pr,ps,rl,rr,res,rs, &
593  rus,rvs,rws,sl,sm,sr,term,ul,ur,us,vl,vml,vmr,vms,vr,vs, &
594  wl,wr,ws,wt1,wt2
595  REAL(RFREAL) :: flx(5)
596  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
597  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,pgv,prhs,psd
598  TYPE(t_global), POINTER :: global
599  TYPE(t_grid), POINTER :: pgrid
600  TYPE(t_patch), POINTER :: ppatch
601 
602 #ifdef SPEC
603  INTEGER :: ispec
604  REAL(RFREAL) :: cps,gcs,mms,y1,y2,ys
605  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcvspec
606  TYPE(t_spec_type), POINTER :: pspectype
607 #endif
608 
609 ! ******************************************************************************
610 ! Start
611 ! ******************************************************************************
612 
613  global => pregion%global
614 
615  CALL registerfunction(global,'RFLU_HLLC_ComputeFlux1_MTCP',&
616  'RFLU_ModHLLCFlux.F90')
617 
618 #ifdef ROCPROF
619  CALL fprofiler_begins("RFLU::HLLC_ComputeFlux1_MTCP")
620 #endif
621 
622 ! ******************************************************************************
623 ! Set dimensions and pointers
624 ! ******************************************************************************
625 
626  pgrid => pregion%grid
627 
628  indgs = pgrid%indGs
629 
630  indcp = pregion%mixtInput%indCp
631  indmf = pregion%mixtInput%indMfMixt
632  indmol = pregion%mixtInput%indMol
633  indsd = pregion%mixtInput%indSd
634 
635  pcv => pregion%mixt%cv
636  pdv => pregion%mixt%dv
637  pgv => pregion%mixt%gv
638 
639  pmf => pregion%mixt%mfMixt
640  prhs => pregion%mixt%rhs
641  psd => pregion%mixt%sd
642 
643 #ifdef SPEC
644  pcvspec => pregion%spec%cv
645 #endif
646 
647  IF ( pregion%spec%cvState /= cv_mixt_state_cons ) THEN
648  CALL errorstop(global,err_cv_state_invalid,__line__)
649  END IF ! pRegion%spec%cvState
650 
651  IF ( pregion%mixtInput%gasModel /= gas_model_mixt_tcperf ) THEN
652  CALL errorstop(global,err_gasmodel_invalid,__line__)
653  END IF ! pRegion%mixtInput%gasModel
654 
655  IF ( (indcp /= 1) .OR. (indmol /= 1) ) THEN
656  CALL errorstop(global,err_reached_default,__line__)
657  END IF ! indCp
658 
659 ! ******************************************************************************
660 ! Compute fluxes through interior faces
661 ! ******************************************************************************
662 
663  DO ifg = 1,pgrid%nFaces
664  c1 = pgrid%f2c(1,ifg)
665  c2 = pgrid%f2c(2,ifg)
666 
667 ! ==============================================================================
668 ! Get face geometry and grid speed
669 ! ==============================================================================
670 
671  nx = pgrid%fn(xcoord,ifg)
672  ny = pgrid%fn(ycoord,ifg)
673  nz = pgrid%fn(zcoord,ifg)
674  nm = pgrid%fn(xyzmag,ifg)
675 
676  fs = pgrid%gs(indgs*ifg)
677 
678 ! ==============================================================================
679 ! Compute left and right states
680 ! ==============================================================================
681 
682 ! ------------------------------------------------------------------------------
683 ! Left state
684 ! ------------------------------------------------------------------------------
685 
686  rl = pcv(cv_mixt_dens,c1)
687  irl = 1.0_rfreal/rl
688 
689  mml = pgv(gv_mixt_mol,c1)
690  cpl = pgv(gv_mixt_cp ,c1)
691 
692  gcl = mixtperf_r_m(mml)
693  gl = mixtperf_g_cpr(cpl,gcl)
694 
695  ul = pcv(cv_mixt_xmom,c1)*irl
696  vl = pcv(cv_mixt_ymom,c1)*irl
697  wl = pcv(cv_mixt_zmom,c1)*irl
698  pl = pdv(dv_mixt_pres,c1)
699  al = pdv(dv_mixt_soun,c1)
700 
701  el = mixtperf_eo_dgpuvw(rl,gl,pl,ul,vl,wl)
702 
703  vml = ul*ul + vl*vl + wl*wl
704  ql = ul*nx + vl*ny + wl*nz - fs
705  hl = el + irl*pl
706 
707 ! ------------------------------------------------------------------------------
708 ! Right state
709 ! ------------------------------------------------------------------------------
710 
711  rr = pcv(cv_mixt_dens,c2)
712  irr = 1.0_rfreal/rr
713 
714  mmr = pgv(gv_mixt_mol,c2)
715  cpr = pgv(gv_mixt_cp ,c2)
716 
717  gcr = mixtperf_r_m(mmr)
718  gr = mixtperf_g_cpr(cpr,gcr)
719 
720  ur = pcv(cv_mixt_xmom,c2)*irr
721  vr = pcv(cv_mixt_ymom,c2)*irr
722  wr = pcv(cv_mixt_zmom,c2)*irr
723  pr = pdv(dv_mixt_pres,c2)
724  ar = pdv(dv_mixt_soun,c2)
725 
726  er = mixtperf_eo_dgpuvw(rr,gr,pr,ur,vr,wr)
727 
728  vmr = ur*ur + vr*vr + wr*wr
729  qr = ur*nx + vr*ny + wr*nz - fs
730  hr = er + irr*pr
731 
732 ! ==============================================================================
733 ! Compute weights
734 ! ==============================================================================
735 
736  wt1 = sqrt(rr/rl)
737  wt2 = 1.0_rfreal/(1.0_rfreal + wt1)
738 
739 ! ==============================================================================
740 ! Compute averages. NOTE this average differs from that of the TCP case in
741 ! that need to recompute gas properties to compute speed of sound. Choose to
742 ! apply Roe-average also to mass fractions, although there is no analysis to
743 ! back this up. It is not clear what is the best way to do this.
744 ! ==============================================================================
745 
746  us = wt2*(ul + wt1*ur)
747  vs = wt2*(vl + wt1*vr)
748  ws = wt2*(wl + wt1*wr)
749  qs = wt2*(ql + wt1*qr)
750  hs = wt2*(hl + wt1*hr)
751 
752 #ifdef SPEC
753  mms = 0.0_rfreal
754  cps = 0.0_rfreal
755 
756  DO ispec = 1,pregion%specInput%nSpecies
757  pspectype => pregion%specInput%specType(ispec)
758 
759  y1 = irl*pcvspec(ispec,c1)
760  y2 = irr*pcvspec(ispec,c2)
761 
762  ys = wt2*(y1 + wt1*y2)
763 
764  mms = mms + ys/pspectype%pMaterial%molw
765  cps = cps + ys*pspectype%pMaterial%spht
766  END DO ! iSpec
767 
768  mms = 1.0_rfreal/mms
769  gcs = mixtperf_r_m(mms)
770  gs = mixtperf_g_cpr(cps,gcs)
771 #endif
772 
773  vms = us*us + vs*vs + ws*ws
774  as = mixtperf_c_ghovm2(gs,hs,vms)
775  sl = min(ql-al,qs-as)
776  sr = max(qr+ar,qs+as)
777 
778 ! ==============================================================================
779 ! Compute fluxes
780 ! ==============================================================================
781 
782  IF ( sl > 0.0_rfreal ) THEN
783  flx(1) = ql* rl *nm
784  flx(2) = (ql* rl*ul + pl*nx )*nm
785  flx(3) = (ql* rl*vl + pl*ny )*nm
786  flx(4) = (ql* rl*wl + pl*nz )*nm
787  flx(5) = (ql*(rl*el + pl) + pl*fs)*nm
788  ELSE IF ( sr < 0.0_rfreal ) THEN
789  flx(1) = qr* rr *nm
790  flx(2) = (qr* rr*ur + pr*nx )*nm
791  flx(3) = (qr* rr*vr + pr*ny )*nm
792  flx(4) = (qr* rr*wr + pr*nz )*nm
793  flx(5) = (qr*(rr*er + pr) + pr*fs)*nm
794  ELSE
795  sm = (rr*qr*(sr-qr) - rl*ql*(sl-ql) - pr + pl)/(rr*(sr-qr) - rl*(sl-ql))
796  ps = pl + rl*(ql-sl)*(ql-sm)
797 
798  IF ( sm > 0.0_rfreal ) THEN
799  term = 1.0_rfreal/(sl-sm)
800 
801  rs = term* (sl-ql)*rl
802  rus = term*((sl-ql)*rl*ul + (ps - pl )*nx)
803  rvs = term*((sl-ql)*rl*vl + (ps - pl )*ny)
804  rws = term*((sl-ql)*rl*wl + (ps - pl )*nz)
805  res = term*((sl-ql)*rl*el + (ps*sm - pl*ql) )
806  ELSE
807  term = 1.0_rfreal/(sr-sm)
808 
809  rs = term* (sr-qr)*rr
810  rus = term*((sr-qr)*rr*ur + (ps - pr )*nx)
811  rvs = term*((sr-qr)*rr*vr + (ps - pr )*ny)
812  rws = term*((sr-qr)*rr*wr + (ps - pr )*nz)
813  res = term*((sr-qr)*rr*er + (ps*sm - pr*qr) )
814  END IF ! sm
815 
816  flx(1) = sm* rs *nm
817  flx(2) = (sm* rus + ps*nx )*nm
818  flx(3) = (sm* rvs + ps*ny )*nm
819  flx(4) = (sm* rws + ps*nz )*nm
820  flx(5) = (sm*(res + ps) - ps*fs)*nm
821  END IF ! sl
822 
823 ! ==============================================================================
824 ! Store mass flux
825 ! ==============================================================================
826 
827  pmf(indmf*ifg) = flx(1)
828 
829 ! ==============================================================================
830 ! Accumulate into residual
831 ! ==============================================================================
832 
833  prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + flx(1)
834  prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + flx(2)
835  prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + flx(3)
836  prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + flx(4)
837  prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + flx(5)
838 
839  prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - flx(1)
840  prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - flx(2)
841  prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - flx(3)
842  prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - flx(4)
843  prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - flx(5)
844 
845 ! ==============================================================================
846 ! Accumulate into substantial derivative
847 ! ==============================================================================
848 
849  psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + us*pmf(indmf*ifg)
850  psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vs*pmf(indmf*ifg)
851  psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + ws*pmf(indmf*ifg)
852 
853  psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - us*pmf(indmf*ifg)
854  psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vs*pmf(indmf*ifg)
855  psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - ws*pmf(indmf*ifg)
856  END DO ! ifg
857 
858 ! ******************************************************************************
859 ! Compute fluxes through boundary faces
860 ! ******************************************************************************
861 
862  DO ipatch = 1,pregion%grid%nPatches
863  ppatch => pregion%patches(ipatch)
864 
865  CALL rflu_centralfirstpatch(pregion,ppatch)
866  END DO ! iPatch
867 
868 ! ******************************************************************************
869 ! End
870 ! ******************************************************************************
871 
872 #ifdef ROCPROF
873  CALL fprofiler_ends("RFLU::HLLC_ComputeFlux1_MTCP")
874 #endif
875 
876  CALL deregisterfunction(global)
877 
878 END SUBROUTINE rflu_hllc_computeflux1_mtcp
879 
880 
881 
882 
883 
884 
885 
886 
887 ! ******************************************************************************
888 !
889 ! Purpose: Compute convective fluxes using first-order accurate version of
890 ! HLLC scheme for thermally and calorically perfect gas.
891 !
892 ! Description: None.
893 !
894 ! Input:
895 ! pRegion Pointer to region
896 !
897 ! Output: None.
898 !
899 ! Notes:
900 ! 1. Only applicable to thermally and calorically perfect gas because would
901 ! need to recompute mixture properties at faces depending on face states.
902 !
903 ! ******************************************************************************
904 
905 SUBROUTINE rflu_hllc_computeflux1_tcp(pRegion)
906 
907  USE modinterfaces, ONLY: mixtperf_c_ghovm2, &
908  rflu_centralfirstpatch
909 
910  IMPLICIT NONE
911 
912 ! ******************************************************************************
913 ! Definitions and declarations
914 ! ******************************************************************************
915 
916 ! ==============================================================================
917 ! Arguments
918 ! ==============================================================================
919 
920  TYPE(t_region), POINTER :: pregion
921 
922 ! ==============================================================================
923 ! Locals
924 ! ==============================================================================
925 
926  INTEGER :: c1,c2,ifg,indgs,indmf,indsd,ipatch
927  INTEGER, DIMENSION(:,:), POINTER :: f2c
928  REAL(RFREAL) :: al,ar,as,el,er,fs,g,hl,hr,hs,irl,irr,nm,nx,ny,nz,ql,qr, &
929  qs,pl,pr,ps,rl,rr,res,rs,rus,rvs,rws,sl,sm,sr,term,ul,ur, &
930  us,vl,vml,vmr,vms,vr,vs,wl,wr,ws,wt1,wt2
931  REAL(RFREAL) :: flx(5)
932  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
933  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,prhs,psd
934  TYPE(t_global), POINTER :: global
935  TYPE(t_grid), POINTER :: pgrid
936  TYPE(t_patch), POINTER :: ppatch
937 
938 ! ******************************************************************************
939 ! Start
940 ! ******************************************************************************
941 
942  global => pregion%global
943 
944  CALL registerfunction(global,'RFLU_HLLC_ComputeFlux1_TCP',&
945  'RFLU_ModHLLCFlux.F90')
946 
947 #ifdef ROCPROF
948  CALL fprofiler_begins("RFLU::HLLC_ComputeFlux1_TCP")
949 #endif
950 
951 ! ******************************************************************************
952 ! Set dimensions and pointers
953 ! ******************************************************************************
954 
955  pgrid => pregion%grid
956 
957  indgs = pgrid%indGs
958 
959  indmf = pregion%mixtInput%indMfMixt
960  indsd = pregion%mixtInput%indSd
961 
962  pcv => pregion%mixt%cv
963  pdv => pregion%mixt%dv
964 
965  pmf => pregion%mixt%mfMixt
966  prhs => pregion%mixt%rhs
967  psd => pregion%mixt%sd
968 
969  IF ( pregion%mixtInput%gasModel /= gas_model_tcperf ) THEN
970  CALL errorstop(global,err_gasmodel_invalid,__line__)
971  END IF ! pRegion%mixtInput%gasModel
972 
973  g = global%refGamma
974 
975 ! ******************************************************************************
976 ! Compute fluxes through interior faces
977 ! ******************************************************************************
978 
979  DO ifg = 1,pgrid%nFaces
980  c1 = pgrid%f2c(1,ifg)
981  c2 = pgrid%f2c(2,ifg)
982 
983 ! ==============================================================================
984 ! Get face geometry and grid speed
985 ! ==============================================================================
986 
987  nx = pgrid%fn(xcoord,ifg)
988  ny = pgrid%fn(ycoord,ifg)
989  nz = pgrid%fn(zcoord,ifg)
990  nm = pgrid%fn(xyzmag,ifg)
991 
992  fs = pgrid%gs(indgs*ifg)
993 
994 ! ==============================================================================
995 ! Compute left and right states
996 ! ==============================================================================
997 
998 ! ------------------------------------------------------------------------------
999 ! Left state
1000 ! ------------------------------------------------------------------------------
1001 
1002  rl = pcv(cv_mixt_dens,c1)
1003  irl = 1.0_rfreal/rl
1004 
1005  ul = pcv(cv_mixt_xmom,c1)*irl
1006  vl = pcv(cv_mixt_ymom,c1)*irl
1007  wl = pcv(cv_mixt_zmom,c1)*irl
1008  el = pcv(cv_mixt_ener,c1)*irl
1009  pl = pdv(dv_mixt_pres,c1)
1010  al = pdv(dv_mixt_soun,c1)
1011 
1012  vml = ul*ul + vl*vl + wl*wl
1013  ql = ul*nx + vl*ny + wl*nz - fs
1014  hl = el + irl*pl
1015 
1016 ! ------------------------------------------------------------------------------
1017 ! Right state
1018 ! ------------------------------------------------------------------------------
1019 
1020  rr = pcv(cv_mixt_dens,c2)
1021  irr = 1.0_rfreal/rr
1022 
1023  ur = pcv(cv_mixt_xmom,c2)*irr
1024  vr = pcv(cv_mixt_ymom,c2)*irr
1025  wr = pcv(cv_mixt_zmom,c2)*irr
1026  er = pcv(cv_mixt_ener,c2)*irr
1027  pr = pdv(dv_mixt_pres,c2)
1028  ar = pdv(dv_mixt_soun,c2)
1029 
1030  vmr = ur*ur + vr*vr + wr*wr
1031  qr = ur*nx + vr*ny + wr*nz - fs
1032  hr = er + irr*pr
1033 
1034 ! ==============================================================================
1035 ! Compute weights
1036 ! ==============================================================================
1037 
1038  wt1 = sqrt(rr/rl)
1039  wt2 = 1.0_rfreal/(1.0_rfreal + wt1)
1040 
1041 ! ==============================================================================
1042 ! Compute averages
1043 ! ==============================================================================
1044 
1045  us = wt2*(ul + wt1*ur)
1046  vs = wt2*(vl + wt1*vr)
1047  ws = wt2*(wl + wt1*wr)
1048  qs = wt2*(ql + wt1*qr)
1049  hs = wt2*(hl + wt1*hr)
1050 
1051  vms = us*us + vs*vs + ws*ws
1052  as = mixtperf_c_ghovm2(g,hs,vms)
1053  sl = min(ql-al,qs-as)
1054  sr = max(qr+ar,qs+as)
1055 
1056 ! ==============================================================================
1057 ! Compute fluxes
1058 ! ==============================================================================
1059 
1060  IF ( sl > 0.0_rfreal ) THEN
1061  flx(1) = ql* rl *nm
1062  flx(2) = (ql* rl*ul + pl*nx )*nm
1063  flx(3) = (ql* rl*vl + pl*ny )*nm
1064  flx(4) = (ql* rl*wl + pl*nz )*nm
1065  flx(5) = (ql*(rl*el + pl) + pl*fs)*nm
1066  ELSE IF ( sr < 0.0_rfreal ) THEN
1067  flx(1) = qr* rr *nm
1068  flx(2) = (qr* rr*ur + pr*nx )*nm
1069  flx(3) = (qr* rr*vr + pr*ny )*nm
1070  flx(4) = (qr* rr*wr + pr*nz )*nm
1071  flx(5) = (qr*(rr*er + pr) + pr*fs)*nm
1072  ELSE
1073  sm = (rr*qr*(sr-qr) - rl*ql*(sl-ql) - pr + pl)/(rr*(sr-qr) - rl*(sl-ql))
1074  ps = pl + rl*(ql-sl)*(ql-sm)
1075 
1076  IF ( sm > 0.0_rfreal ) THEN
1077  term = 1.0_rfreal/(sl-sm)
1078 
1079  rs = term* (sl-ql)*rl
1080  rus = term*((sl-ql)*rl*ul + (ps - pl )*nx)
1081  rvs = term*((sl-ql)*rl*vl + (ps - pl )*ny)
1082  rws = term*((sl-ql)*rl*wl + (ps - pl )*nz)
1083  res = term*((sl-ql)*rl*el + (ps*sm - pl*ql) )
1084  ELSE
1085  term = 1.0_rfreal/(sr-sm)
1086 
1087  rs = term* (sr-qr)*rr
1088  rus = term*((sr-qr)*rr*ur + (ps - pr )*nx)
1089  rvs = term*((sr-qr)*rr*vr + (ps - pr )*ny)
1090  rws = term*((sr-qr)*rr*wr + (ps - pr )*nz)
1091  res = term*((sr-qr)*rr*er + (ps*sm - pr*qr) )
1092  END IF ! sm
1093 
1094  flx(1) = sm* rs *nm
1095  flx(2) = (sm* rus + ps*nx )*nm
1096  flx(3) = (sm* rvs + ps*ny )*nm
1097  flx(4) = (sm* rws + ps*nz )*nm
1098  flx(5) = (sm*(res + ps) - ps*fs)*nm
1099  END IF ! sl
1100 
1101 ! ==============================================================================
1102 ! Store mass flux
1103 ! ==============================================================================
1104 
1105  pmf(indmf*ifg) = flx(1)
1106 
1107 ! ==============================================================================
1108 ! Accumulate into residual
1109 ! ==============================================================================
1110 
1111  prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + flx(1)
1112  prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + flx(2)
1113  prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + flx(3)
1114  prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + flx(4)
1115  prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + flx(5)
1116 
1117  prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - flx(1)
1118  prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - flx(2)
1119  prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - flx(3)
1120  prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - flx(4)
1121  prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - flx(5)
1122 
1123 ! ==============================================================================
1124 ! Accumulate into substantial derivative
1125 ! ==============================================================================
1126 
1127  psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + us*pmf(indmf*ifg)
1128  psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vs*pmf(indmf*ifg)
1129  psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + ws*pmf(indmf*ifg)
1130 
1131  psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - us*pmf(indmf*ifg)
1132  psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vs*pmf(indmf*ifg)
1133  psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - ws*pmf(indmf*ifg)
1134  END DO ! ifg
1135 
1136 ! ******************************************************************************
1137 ! Compute fluxes through boundary faces
1138 ! ******************************************************************************
1139 
1140  DO ipatch = 1,pregion%grid%nPatches
1141  ppatch => pregion%patches(ipatch)
1142 
1143  CALL rflu_centralfirstpatch(pregion,ppatch)
1144  END DO ! iPatch
1145 
1146 ! ******************************************************************************
1147 ! End
1148 ! ******************************************************************************
1149 
1150 #ifdef ROCPROF
1151  CALL fprofiler_ends("RFLU::HLLC_ComputeFlux1_TCP")
1152 #endif
1153 
1154  CALL deregisterfunction(global)
1155 
1156 END SUBROUTINE rflu_hllc_computeflux1_tcp
1157 
1158 
1159 
1160 
1161 
1162 
1163 
1164 
1165 ! ******************************************************************************
1166 !
1167 ! Purpose: Compute central convective fluxes for multiphase mixture using
1168 ! second-order accurate version of HLLC scheme.
1169 !
1170 ! Description: None.
1171 !
1172 ! Input:
1173 ! pRegion Data of current region.
1174 !
1175 ! Output: None.
1176 !
1177 ! Notes: None.
1178 !
1179 ! ******************************************************************************
1180 
1181 SUBROUTINE rflu_hllc_computeflux2_gl(pRegion)
1182 
1183  USE modinterfaces, ONLY: mixtgasliq_c, &
1184  mixtgasliq_p, &
1185  mixtliq_d_dobpppobttto, &
1186  mixtliq_c2_bp, &
1187  mixtperf_c2_grt, &
1188  mixtperf_cv_cpr, &
1189  mixtperf_d_prt, &
1190  mixtperf_r_cpg, &
1191  mixtperf_r_m, &
1192  rflu_centralfirstpatch_gl, &
1193  rflu_centralsecondpatch_gl
1194 
1195  IMPLICIT NONE
1196 
1197 ! ******************************************************************************
1198 ! Definitions and declarations
1199 ! ******************************************************************************
1200 
1201 ! ==============================================================================
1202 ! Arguments
1203 ! ==============================================================================
1204 
1205  TYPE(t_region), POINTER :: pregion
1206 
1207 ! ==============================================================================
1208 ! Locals
1209 ! ==============================================================================
1210 
1211  INTEGER :: c1,c2,ifg,indgs,indmf,indsd,ipatch
1212  REAL(RFREAL) :: bgh2,bgl2,bgr2,blh2,bll2,blr2,bp,bt,bvh2,bvl2,bvr2,cml, &
1213  cmr,cms,cvg,cvl,cvm,cvv,cgh2,cgl2,cgr2,clh2,cll2,clr2, &
1214  cvh2,cvl2,cvr2,dx,dy,dz,el,er,fs,irl,irr,nm,nx,ny,nz,ph, &
1215  pl,po,pr,ps,ql,qr,qs,rel,rer,res,rgh,rgl,rgr,rh,rl,rlh, &
1216  rll,rlr,ro,rr,rs,rus,rvh,rvl,rvr,rvs,rws,rygl,rygr,ryll, &
1217  rylr,ryvl,ryvr,rg,rv,sl,sm,sr,term,th,tl,to,tr,ul,ur,us, &
1218  vfgh,vfgl,vfgr,vflh,vfll,vflr,vfvh,vfvl,vfvr,vl,vr,vs,vl2, &
1219  vms,vr2,wl,wr,ws,wt1,wt2,xc,yc,zc,ygl,ygr,yll,ylr,yvl,yvr
1220  REAL(RFREAL) :: flx(5)
1221  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
1222  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcvmixt,pcvspec,pdvmixt,prhs,psd
1223  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgradmixt,pgradspec
1224  TYPE(t_global), POINTER :: global
1225  TYPE(t_grid), POINTER :: pgrid
1226  TYPE(t_patch), POINTER :: ppatch
1227 
1228 ! ******************************************************************************
1229 ! Start
1230 ! ******************************************************************************
1231 
1232  global => pregion%global
1233 
1234  CALL registerfunction(global,'RFLU_HLLC_ComputeFlux2_GL',&
1235  'RFLU_ModHLLCFlux.F90')
1236 
1237 ! ******************************************************************************
1238 ! Set dimensions and pointers
1239 ! ******************************************************************************
1240 
1241  pgrid => pregion%grid
1242 
1243  indgs = pgrid%indGs
1244 
1245  indmf = pregion%mixtInput%indMfMixt
1246  indsd = pregion%mixtInput%indSd
1247 
1248  pcvmixt => pregion%mixt%cv
1249  pdvmixt => pregion%mixt%dv
1250  pcvspec => pregion%spec%cv
1251 
1252  pgradmixt => pregion%mixt%gradCell
1253  pgradspec => pregion%spec%gradCell
1254 
1255  pmf => pregion%mixt%mfMixt
1256  prhs => pregion%mixt%rhs
1257  psd => pregion%mixt%sd
1258 
1259  IF ( pregion%spec%cvState /= cv_mixt_state_cons ) THEN
1260  CALL errorstop(global,err_cv_state_invalid,__line__)
1261  END IF ! pRegion%spec%cvState
1262 
1263  IF ( pregion%mixtInput%gasModel /= gas_model_mixt_gasliq ) THEN
1264  CALL errorstop(global,err_gasmodel_invalid,__line__)
1265  END IF ! pRegion%mixtInput%gasModel
1266 
1267 ! ******************************************************************************
1268 ! Define constants
1269 ! ******************************************************************************
1270 
1271  ro = global%refDensityLiq
1272  po = global%refPressLiq
1273  to = global%refTempLiq
1274  bp = global%refBetaPLiq
1275  bt = global%refBetaTLiq
1276  cvl = global%refCvLiq
1277 
1278  rg = mixtperf_r_m(pregion%specInput%specType(1)%pMaterial%molw)
1279  cvg = mixtperf_cv_cpr(pregion%specInput%specType(1)%pMaterial%spht,rg)
1280 
1281  rv = mixtperf_r_m(pregion%specInput%specType(2)%pMaterial%molw)
1282  cvv = mixtperf_cv_cpr(pregion%specInput%specType(2)%pMaterial%spht,rv)
1283 
1284 ! *****************************************************************************
1285 ! Compute fluxes through interior faces
1286 ! *****************************************************************************
1287 
1288  DO ifg = 1,pgrid%nFaces
1289  c1 = pgrid%f2c(1,ifg)
1290  c2 = pgrid%f2c(2,ifg)
1291 
1292 ! =============================================================================
1293 ! Get face geometry and grid speed
1294 ! =============================================================================
1295 
1296  nx = pgrid%fn(xcoord,ifg)
1297  ny = pgrid%fn(ycoord,ifg)
1298  nz = pgrid%fn(zcoord,ifg)
1299  nm = pgrid%fn(xyzmag,ifg)
1300 
1301  xc = pgrid%fc(xcoord,ifg)
1302  yc = pgrid%fc(ycoord,ifg)
1303  zc = pgrid%fc(zcoord,ifg)
1304 
1305  fs = pgrid%gs(indgs*ifg)
1306 
1307 ! =============================================================================
1308 ! Compute left and right states
1309 ! =============================================================================
1310 
1311 ! -----------------------------------------------------------------------------
1312 ! Left state
1313 ! -----------------------------------------------------------------------------
1314 
1315  rl = pcvmixt(cv_mixt_dens,c1)
1316  irl = 1.0_rfreal/rl
1317 
1318  ul = pcvmixt(cv_mixt_xmom,c1)*irl
1319  vl = pcvmixt(cv_mixt_ymom,c1)*irl
1320  wl = pcvmixt(cv_mixt_zmom,c1)*irl
1321  tl = pdvmixt(dv_mixt_temp,c1)
1322 
1323  ygl = pcvspec(1,c1)*irl
1324  yvl = pcvspec(2,c1)*irl
1325 
1326  dx = xc - pgrid%cofg(xcoord,c1)
1327  dy = yc - pgrid%cofg(ycoord,c1)
1328  dz = zc - pgrid%cofg(zcoord,c1)
1329 
1330  rl = rl + pgradmixt(xcoord,grc_mixt_dens,c1)*dx &
1331  + pgradmixt(ycoord,grc_mixt_dens,c1)*dy &
1332  + pgradmixt(zcoord,grc_mixt_dens,c1)*dz
1333  ul = ul + pgradmixt(xcoord,grc_mixt_xvel,c1)*dx &
1334  + pgradmixt(ycoord,grc_mixt_xvel,c1)*dy &
1335  + pgradmixt(zcoord,grc_mixt_xvel,c1)*dz
1336  vl = vl + pgradmixt(xcoord,grc_mixt_yvel,c1)*dx &
1337  + pgradmixt(ycoord,grc_mixt_yvel,c1)*dy &
1338  + pgradmixt(zcoord,grc_mixt_yvel,c1)*dz
1339  wl = wl + pgradmixt(xcoord,grc_mixt_zvel,c1)*dx &
1340  + pgradmixt(ycoord,grc_mixt_zvel,c1)*dy &
1341  + pgradmixt(zcoord,grc_mixt_zvel,c1)*dz
1342  tl = tl + pgradmixt(xcoord,grc_mixt_temp,c1)*dx &
1343  + pgradmixt(ycoord,grc_mixt_temp,c1)*dy &
1344  + pgradmixt(zcoord,grc_mixt_temp,c1)*dz
1345 
1346  ygl = ygl + pgradspec(xcoord,1,c1)*dx &
1347  + pgradspec(ycoord,1,c1)*dy &
1348  + pgradspec(zcoord,1,c1)*dz
1349  yvl = yvl + pgradspec(xcoord,2,c1)*dx &
1350  + pgradspec(ycoord,2,c1)*dy &
1351  + pgradspec(zcoord,2,c1)*dz
1352 
1353  IF ( ygl > 1.0_rfreal ) THEN
1354  ygl = 1.0_rfreal
1355  yvl = 0.0_rfreal
1356  ELSE IF ( ygl < 0.0_rfreal ) THEN
1357  ygl = 0.0_rfreal
1358  ELSE IF ( yvl > 1.0_rfreal ) THEN
1359  yvl = 1.0_rfreal
1360  ygl = 0.0_rfreal
1361  ELSE IF ( yvl < 0.0_rfreal ) THEN
1362  yvl = 0.0_rfreal
1363  END IF ! Ygl
1364 
1365  yll = 1.0_rfreal - ygl - yvl
1366 
1367  IF ( yll > 1.0_rfreal) THEN
1368  yll = 1.0_rfreal
1369  ELSE IF ( yll < 0.0_rfreal) THEN
1370  yll = 0.0_rfreal
1371  END IF ! Yll
1372 
1373  vl2 = ul*ul + vl*vl + wl*wl
1374  ql = ul*nx + vl*ny + wl*nz - fs
1375 
1376  cll2 = mixtliq_c2_bp(bp)
1377  cvl2 = mixtperf_c2_grt(1.0_rfreal,rv,tl)
1378  cgl2 = mixtperf_c2_grt(1.0_rfreal,rg,tl)
1379 
1380  rygl = rl*ygl
1381  ryvl = rl*yvl
1382  ryll = rl*yll
1383 
1384  pl = mixtgasliq_p(ryll,ryvl,rygl,cll2,cvl2,cgl2,rl,ro,po,to,bp,bt,tl)
1385 
1386  rll = mixtliq_d_dobpppobttto(ro,bp,bt,pl,po,tl,to)
1387  rvl = mixtperf_d_prt(pl,rv,tl)
1388  rgl = mixtperf_d_prt(pl,rg,tl)
1389 
1390  vfgl = rygl/rgl
1391  vfvl = ryvl/rvl
1392  vfll = ryll/rll
1393 
1394  cvm = (ryll*cvl + ryvl*cvv + rygl*cvg)/rl
1395  rel = rl*cvm*tl + 0.5*rl*vl2
1396  el = rel/rl
1397 
1398  bll2 = -bt/bp
1399  bvl2 = rvl*rv
1400  bgl2 = rgl*rg
1401 
1402  cml = mixtgasliq_c(cvm,rl,pl,rll,rvl,rgl,vfll,vfvl,vfgl,cll2,cvl2,cgl2, &
1403  bll2,bvl2,bgl2)
1404 
1405 ! -----------------------------------------------------------------------------
1406 ! Right state
1407 ! -----------------------------------------------------------------------------
1408 
1409  rr = pcvmixt(cv_mixt_dens,c2)
1410  irr = 1.0_rfreal/rr
1411 
1412  ur = pcvmixt(cv_mixt_xmom,c2)*irr
1413  vr = pcvmixt(cv_mixt_ymom,c2)*irr
1414  wr = pcvmixt(cv_mixt_zmom,c2)*irr
1415  tr = pdvmixt(dv_mixt_temp,c2)
1416 
1417  ygr = pcvspec(1,c2)*irr
1418  yvr = pcvspec(2,c2)*irr
1419 
1420  dx = xc - pgrid%cofg(xcoord,c2)
1421  dy = yc - pgrid%cofg(ycoord,c2)
1422  dz = zc - pgrid%cofg(zcoord,c2)
1423 
1424  rr = rr + pgradmixt(xcoord,grc_mixt_dens,c2)*dx &
1425  + pgradmixt(ycoord,grc_mixt_dens,c2)*dy &
1426  + pgradmixt(zcoord,grc_mixt_dens,c2)*dz
1427  ur = ur + pgradmixt(xcoord,grc_mixt_xvel,c2)*dx &
1428  + pgradmixt(ycoord,grc_mixt_xvel,c2)*dy &
1429  + pgradmixt(zcoord,grc_mixt_xvel,c2)*dz
1430  vr = vr + pgradmixt(xcoord,grc_mixt_yvel,c2)*dx &
1431  + pgradmixt(ycoord,grc_mixt_yvel,c2)*dy &
1432  + pgradmixt(zcoord,grc_mixt_yvel,c2)*dz
1433  wr = wr + pgradmixt(xcoord,grc_mixt_zvel,c2)*dx &
1434  + pgradmixt(ycoord,grc_mixt_zvel,c2)*dy &
1435  + pgradmixt(zcoord,grc_mixt_zvel,c2)*dz
1436  tr = tr + pgradmixt(xcoord,grc_mixt_temp,c2)*dx &
1437  + pgradmixt(ycoord,grc_mixt_temp,c2)*dy &
1438  + pgradmixt(zcoord,grc_mixt_temp,c2)*dz
1439 
1440  ygr = ygr + pgradspec(xcoord,1,c2)*dx &
1441  + pgradspec(ycoord,1,c2)*dy &
1442  + pgradspec(zcoord,1,c2)*dz
1443  yvr = yvr + pgradspec(xcoord,2,c2)*dx &
1444  + pgradspec(ycoord,2,c2)*dy &
1445  + pgradspec(zcoord,2,c2)*dz
1446 
1447  IF ( ygr > 1.0_rfreal ) THEN
1448  ygr = 1.0_rfreal
1449  yvr = 0.0_rfreal
1450  ELSE IF ( ygr < 0.0_rfreal ) THEN
1451  ygr = 0.0_rfreal
1452  ELSE IF ( yvr > 1.0_rfreal ) THEN
1453  yvr = 1.0_rfreal
1454  ygr = 0.0_rfreal
1455  ELSE IF ( yvr < 0.0_rfreal ) THEN
1456  yvr = 0.0_rfreal
1457  END IF ! Ygr
1458 
1459  ylr = 1.0_rfreal - ygr - yvr
1460 
1461  IF ( ylr > 1.0_rfreal ) THEN
1462  ylr = 1.0_rfreal
1463  ELSE IF ( ylr < 0.0_rfreal ) THEN
1464  ylr = 0.0_rfreal
1465  END IF ! Ylr
1466 
1467  vr2 = ur*ur + vr*vr + wr*wr
1468  qr = ur*nx + vr*ny + wr*nz - fs
1469 
1470  clr2 = mixtliq_c2_bp(bp)
1471  cvr2 = mixtperf_c2_grt(1.0_rfreal,rv,tr)
1472  cgr2 = mixtperf_c2_grt(1.0_rfreal,rg,tr)
1473 
1474  rygr = rr*ygr
1475  ryvr = rr*yvr
1476  rylr = rr*ylr
1477 
1478  pr = mixtgasliq_p(rylr,ryvr,rygr,clr2,cvr2,cgr2,rr,ro,po,to,bp,bt,tr)
1479 
1480  rlr = mixtliq_d_dobpppobttto(ro,bp,bt,pr,po,tr,to)
1481  rvr = mixtperf_d_prt(pr,rv,tr)
1482  rgr = mixtperf_d_prt(pr,rg,tr)
1483 
1484  vfgr = rygr/rgr
1485  vfvr = ryvr/rvr
1486  vflr = rylr/rlr
1487 
1488  cvm = (rylr*cvl + ryvr*cvv + rygr*cvg)/rr
1489  rer = rr*cvm*tr + 0.5*rr*vr2
1490  er = rer/rr
1491 
1492  blr2 = -bt/bp
1493  bvr2 = rvr*rv
1494  bgr2 = rgr*rg
1495 
1496  cmr = mixtgasliq_c(cvm,rr,pr,rlr,rvr,rgr,vflr,vfvr,vfgr,clr2,cvr2,cgr2, &
1497  blr2,bvr2,bgr2)
1498 
1499 ! ==============================================================================
1500 ! Compute weights
1501 ! ==============================================================================
1502 
1503  wt1 = sqrt(rr/rl)
1504  wt2 = 1.0_rfreal/(1.0_rfreal + wt1)
1505 
1506 ! ==============================================================================
1507 ! Compute averages
1508 ! ==============================================================================
1509 
1510  us = wt2*(ul + wt1*ur)
1511  vs = wt2*(vl + wt1*vr)
1512  ws = wt2*(wl + wt1*wr)
1513  qs = wt2*(ql + wt1*qr)
1514 
1515  vms = us*us + vs*vs + ws*ws
1516 
1517  rh = wt1*rl
1518  ph = wt2*(pl + wt1*pr)
1519  th = wt2*(tl + wt1*tr)
1520 
1521  vfgh = wt2*(vfgl + wt1*vfgr)
1522  vfvh = wt2*(vfvl + wt1*vfvr)
1523 
1524  IF ( vfgh > 1.0_rfreal ) THEN
1525  vfgh = 1.0_rfreal
1526  vfvh = 0.0_rfreal
1527  ELSE IF ( vfgh < 0.0_rfreal ) THEN
1528  vfgh = 0.0_rfreal
1529  ELSE IF ( vfvh > 1.0_rfreal ) THEN
1530  vfvh = 1.0_rfreal
1531  vfgh = 0.0_rfreal
1532  ELSE IF ( vfvh < 0.0_rfreal ) THEN
1533  vfvh = 0.0_rfreal
1534  END IF ! vgfh
1535 
1536  vflh = 1.0_rfreal - vfvh - vfgh
1537 
1538  IF ( vflh > 1.0_rfreal ) THEN
1539  vflh = 1.0_rfreal
1540  ELSE IF ( vflh < 0.0_rfreal ) THEN
1541  vflh = 0.0_rfreal
1542  END IF ! vflh
1543 
1544  clh2 = mixtliq_c2_bp(bp)
1545  cvh2 = mixtperf_c2_grt(1.0_rfreal,rv,th)
1546  cgh2 = mixtperf_c2_grt(1.0_rfreal,rg,th)
1547 
1548  rlh = mixtliq_d_dobpppobttto(ro,bp,bt,ph,po,th,to)
1549  rvh = mixtperf_d_prt(ph,rv,th)
1550  rgh = mixtperf_d_prt(ph,rg,th)
1551 
1552  blh2 = -bt/bp
1553  bvh2 = rvh*rv
1554  bgh2 = rgh*rg
1555 
1556  cvm = (rlh*vflh*cvl + rvh*vfvh*cvv + rgh*vfgh*cvg)/rh
1557 
1558  cms = mixtgasliq_c(cvm,rh,ph,rlh,rvh,rgh,vflh,vfvh,vfgh,clh2,cvh2,cgh2, &
1559  blh2,bvh2,bgh2)
1560 
1561  sl = min(ql-cml,qs-cms)
1562  sr = max(qr+cmr,qs+cms)
1563 
1564 ! ==============================================================================
1565 ! Compute fluxes
1566 ! ==============================================================================
1567 
1568  IF ( sl > 0.0_rfreal ) THEN
1569  flx(1) = ql* rl *nm
1570  flx(2) = (ql* rl*ul + pl*nx )*nm
1571  flx(3) = (ql* rl*vl + pl*ny )*nm
1572  flx(4) = (ql* rl*wl + pl*nz )*nm
1573  flx(5) = (ql*(rl*el + pl) + pl*fs)*nm
1574  ELSE IF ( sr < 0.0_rfreal ) THEN
1575  flx(1) = qr* rr *nm
1576  flx(2) = (qr* rr*ur + pr*nx )*nm
1577  flx(3) = (qr* rr*vr + pr*ny )*nm
1578  flx(4) = (qr* rr*wr + pr*nz )*nm
1579  flx(5) = (qr*(rr*er + pr) + pr*fs)*nm
1580  ELSE
1581  sm = (rr*qr*(sr-qr) - rl*ql*(sl-ql) - pr + pl)/(rr*(sr-qr) - rl*(sl-ql))
1582  ps = pl + rl*(ql-sl)*(ql-sm)
1583 
1584  IF ( sm > 0.0_rfreal ) THEN
1585  term = 1.0_rfreal/(sl-sm)
1586 
1587  rs = term* (sl-ql)*rl
1588  rus = term*((sl-ql)*rl*ul + (ps - pl )*nx)
1589  rvs = term*((sl-ql)*rl*vl + (ps - pl )*ny)
1590  rws = term*((sl-ql)*rl*wl + (ps - pl )*nz)
1591  res = term*((sl-ql)*rl*el + (ps*sm - pl*ql) )
1592  ELSE
1593  term = 1.0_rfreal/(sr-sm)
1594 
1595  rs = term* (sr-qr)*rr
1596  rus = term*((sr-qr)*rr*ur + (ps - pr )*nx)
1597  rvs = term*((sr-qr)*rr*vr + (ps - pr )*ny)
1598  rws = term*((sr-qr)*rr*wr + (ps - pr )*nz)
1599  res = term*((sr-qr)*rr*er + (ps*sm - pr*qr) )
1600  END IF ! sm
1601 
1602  flx(1) = sm* rs *nm
1603  flx(2) = (sm* rus + ps*nx )*nm
1604  flx(3) = (sm* rvs + ps*ny )*nm
1605  flx(4) = (sm* rws + ps*nz )*nm
1606  flx(5) = (sm*(res + ps) - ps*fs)*nm
1607  END IF ! sl
1608 
1609 ! ==============================================================================
1610 ! Store mass flux
1611 ! ==============================================================================
1612 
1613  pmf(indmf*ifg) = flx(1)
1614 
1615 ! ==============================================================================
1616 ! Accumulate into residual
1617 ! ==============================================================================
1618 
1619  prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + flx(1)
1620  prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + flx(2)
1621  prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + flx(3)
1622  prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + flx(4)
1623  prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + flx(5)
1624 
1625  prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - flx(1)
1626  prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - flx(2)
1627  prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - flx(3)
1628  prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - flx(4)
1629  prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - flx(5)
1630 
1631 ! ==============================================================================
1632 ! Accumulate into substantial derivative
1633 ! ==============================================================================
1634 
1635  psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + us*pmf(indmf*ifg)
1636  psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vs*pmf(indmf*ifg)
1637  psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + ws*pmf(indmf*ifg)
1638 
1639  psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - us*pmf(indmf*ifg)
1640  psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vs*pmf(indmf*ifg)
1641  psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - ws*pmf(indmf*ifg)
1642  END DO ! ifg
1643 
1644 ! ******************************************************************************
1645 ! Compute fluxes through boundary faces
1646 ! ******************************************************************************
1647 
1648  DO ipatch = 1,pregion%grid%nPatches
1649  ppatch => pregion%patches(ipatch)
1650 
1651  SELECT CASE ( ppatch%spaceOrder )
1652  CASE ( 1 )
1653  CALL rflu_centralfirstpatch_gl(pregion,ppatch)
1654  CASE ( 2 )
1655  CALL rflu_centralsecondpatch_gl(pregion,ppatch)
1656  CASE default
1657  CALL errorstop(global,err_reached_default,__line__)
1658  END SELECT ! pPatch%spaceOrder
1659  END DO ! iPatch
1660 
1661 ! ******************************************************************************
1662 ! End
1663 ! ******************************************************************************
1664 
1665  CALL deregisterfunction(global)
1666 
1667 END SUBROUTINE rflu_hllc_computeflux2_gl
1668 
1669 
1670 
1671 
1672 
1673 
1674 
1675 
1676 ! ******************************************************************************
1677 !
1678 ! Purpose: Compute convective fluxes using second-order accurate version of
1679 ! HLLC scheme for mixture of thermally and calorically perfect gases.
1680 !
1681 ! Description: None.
1682 !
1683 ! Input:
1684 ! pRegion Pointer to region
1685 !
1686 ! Output: None.
1687 !
1688 ! Notes: None.
1689 !
1690 ! ******************************************************************************
1691 
1692 SUBROUTINE rflu_hllc_computeflux2_mtcp(pRegion)
1693 
1694 #ifdef SPEC
1695  USE modspecies, ONLY: t_spec_type
1696 #endif
1697 
1698  USE modinterfaces, ONLY: mixtperf_c_dgp, &
1699  mixtperf_c_ghovm2, &
1700  mixtperf_eo_dgpuvw, &
1701  mixtperf_g_cpr, &
1702  mixtperf_r_m, &
1703  rflu_centralfirstpatch, &
1704  rflu_centralsecondpatch
1705 
1706  IMPLICIT NONE
1707 
1708 ! ******************************************************************************
1709 ! Definitions and declarations
1710 ! ******************************************************************************
1711 
1712 ! ==============================================================================
1713 ! Arguments
1714 ! ==============================================================================
1715 
1716  TYPE(t_region), POINTER :: pregion
1717 
1718 ! ==============================================================================
1719 ! Locals
1720 ! ==============================================================================
1721 
1722  INTEGER :: c1,c2,ifg,indcp,indgs,indmf,indmol,indsd,ipatch
1723  REAL(RFREAL) :: al,ar,as,dx1,dx2,dy1,dy2,dz1,dz2,el,er,fs,gcl,gcr,gl,gr,gs, &
1724  hl,hr,hs,irl,irr,nm,nx,ny,nz,ql,qr,qs,pl,pr,ps,rl,rr,res, &
1725  rs,rus,rvs,rws,sl,sm,sr,term,ul,ur,us,vl,vml,vmr,vms,vr,vs, &
1726  wl,wr,ws,wt1,wt2,xc,yc,zc
1727  REAL(RFREAL) :: flx(5)
1728  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
1729  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,pgv,prhs,psd
1730  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgc
1731  TYPE(t_global), POINTER :: global
1732  TYPE(t_grid), POINTER :: pgrid
1733  TYPE(t_patch), POINTER :: ppatch
1734 
1735 #ifdef SPEC
1736  INTEGER :: ispec
1737  REAL(RFREAL) :: cpl,cpr,cps,gcs,mml,mmr,mms,y1,y2,ys
1738  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcvspec
1739  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgcspec
1740  TYPE(t_spec_type), POINTER :: pspectype
1741 #endif
1742 
1743 ! ******************************************************************************
1744 ! Start
1745 ! ******************************************************************************
1746 
1747  global => pregion%global
1748 
1749  CALL registerfunction(global,'RFLU_HLLC_ComputeFlux2_MTCP',&
1750  'RFLU_ModHLLCFlux.F90')
1751 
1752 #ifdef ROCPROF
1753  CALL fprofiler_begins("RFLU::HLLC_ComputeFlux2_MTCP")
1754 #endif
1755 
1756 ! ******************************************************************************
1757 ! Set dimensions and pointers
1758 ! ******************************************************************************
1759 
1760  pgrid => pregion%grid
1761 
1762  indgs = pgrid%indGs
1763 
1764  indcp = pregion%mixtInput%indCp
1765  indmf = pregion%mixtInput%indMfMixt
1766  indmol = pregion%mixtInput%indMol
1767  indsd = pregion%mixtInput%indSd
1768 
1769  pcv => pregion%mixt%cv
1770  pdv => pregion%mixt%dv
1771  pgv => pregion%mixt%gv
1772 
1773  pgc => pregion%mixt%gradCell
1774  pmf => pregion%mixt%mfMixt
1775  prhs => pregion%mixt%rhs
1776  psd => pregion%mixt%sd
1777 
1778 #ifdef SPEC
1779  pcvspec => pregion%spec%cv
1780  pgcspec => pregion%spec%gradCell
1781 #endif
1782 
1783  IF ( pregion%spec%cvState /= cv_mixt_state_cons ) THEN
1784  CALL errorstop(global,err_cv_state_invalid,__line__)
1785  END IF ! pRegion%spec%cvState
1786 
1787  IF ( pregion%mixtInput%gasModel /= gas_model_mixt_tcperf ) THEN
1788  CALL errorstop(global,err_gasmodel_invalid,__line__)
1789  END IF ! pRegion%mixtInput%gasModel
1790 
1791  IF ( (indcp /= 1) .OR. (indmol /= 1) ) THEN
1792  CALL errorstop(global,err_reached_default,__line__)
1793  END IF ! indCp
1794 
1795 ! ******************************************************************************
1796 ! Compute fluxes through interior faces
1797 ! ******************************************************************************
1798 
1799  DO ifg = 1,pgrid%nFaces
1800  c1 = pgrid%f2c(1,ifg)
1801  c2 = pgrid%f2c(2,ifg)
1802 
1803 ! ==============================================================================
1804 ! Get face geometry and grid speed
1805 ! ==============================================================================
1806 
1807  nx = pgrid%fn(xcoord,ifg)
1808  ny = pgrid%fn(ycoord,ifg)
1809  nz = pgrid%fn(zcoord,ifg)
1810  nm = pgrid%fn(xyzmag,ifg)
1811 
1812  xc = pgrid%fc(xcoord,ifg)
1813  yc = pgrid%fc(ycoord,ifg)
1814  zc = pgrid%fc(zcoord,ifg)
1815 
1816  fs = pgrid%gs(indgs*ifg)
1817 
1818 ! ==============================================================================
1819 ! Compute left and right states
1820 ! ==============================================================================
1821 
1822 ! ------------------------------------------------------------------------------
1823 ! Left state
1824 ! ------------------------------------------------------------------------------
1825 
1826  rl = pcv(cv_mixt_dens,c1)
1827  irl = 1.0_rfreal/rl
1828 
1829  dx1 = xc - pgrid%cofg(xcoord,c1)
1830  dy1 = yc - pgrid%cofg(ycoord,c1)
1831  dz1 = zc - pgrid%cofg(zcoord,c1)
1832 
1833 #ifdef SPEC
1834  mml = 0.0_rfreal
1835  cpl = 0.0_rfreal
1836 
1837  DO ispec = 1,pregion%specInput%nSpecies
1838  pspectype => pregion%specInput%specType(ispec)
1839 
1840  y1 = irl*pcvspec(ispec,c1) + pgcspec(xcoord,ispec,c1)*dx1 &
1841  + pgcspec(ycoord,ispec,c1)*dy1 &
1842  + pgcspec(zcoord,ispec,c1)*dz1
1843 
1844  mml = mml + y1/pspectype%pMaterial%molw
1845  cpl = cpl + y1*pspectype%pMaterial%spht
1846  END DO ! iSpec
1847 
1848  mml = 1.0_rfreal/mml
1849  gcl = mixtperf_r_m(mml)
1850  gl = mixtperf_g_cpr(cpl,gcl)
1851 #endif
1852 
1853  ul = pcv(cv_mixt_xmom,c1)*irl
1854  vl = pcv(cv_mixt_ymom,c1)*irl
1855  wl = pcv(cv_mixt_zmom,c1)*irl
1856  pl = pdv(dv_mixt_pres,c1)
1857 
1858  rl = rl + pgc(xcoord,grc_mixt_dens,c1)*dx1 &
1859  + pgc(ycoord,grc_mixt_dens,c1)*dy1 &
1860  + pgc(zcoord,grc_mixt_dens,c1)*dz1
1861  ul = ul + pgc(xcoord,grc_mixt_xvel,c1)*dx1 &
1862  + pgc(ycoord,grc_mixt_xvel,c1)*dy1 &
1863  + pgc(zcoord,grc_mixt_xvel,c1)*dz1
1864  vl = vl + pgc(xcoord,grc_mixt_yvel,c1)*dx1 &
1865  + pgc(ycoord,grc_mixt_yvel,c1)*dy1 &
1866  + pgc(zcoord,grc_mixt_yvel,c1)*dz1
1867  wl = wl + pgc(xcoord,grc_mixt_zvel,c1)*dx1 &
1868  + pgc(ycoord,grc_mixt_zvel,c1)*dy1 &
1869  + pgc(zcoord,grc_mixt_zvel,c1)*dz1
1870  pl = pl + pgc(xcoord,grc_mixt_pres,c1)*dx1 &
1871  + pgc(ycoord,grc_mixt_pres,c1)*dy1 &
1872  + pgc(zcoord,grc_mixt_pres,c1)*dz1
1873 
1874  el = mixtperf_eo_dgpuvw(rl,gl,pl,ul,vl,wl)
1875  al = mixtperf_c_dgp(rl,gl,pl)
1876 
1877  vml = ul*ul + vl*vl + wl*wl
1878  ql = ul*nx + vl*ny + wl*nz - fs
1879  hl = el + irl*pl
1880 
1881 ! ------------------------------------------------------------------------------
1882 ! Right state
1883 ! ------------------------------------------------------------------------------
1884 
1885  rr = pcv(cv_mixt_dens,c2)
1886  irr = 1.0_rfreal/rr
1887 
1888  dx2 = xc - pgrid%cofg(xcoord,c2)
1889  dy2 = yc - pgrid%cofg(ycoord,c2)
1890  dz2 = zc - pgrid%cofg(zcoord,c2)
1891 
1892 #ifdef SPEC
1893  mmr = 0.0_rfreal
1894  cpr = 0.0_rfreal
1895 
1896  DO ispec = 1,pregion%specInput%nSpecies
1897  pspectype => pregion%specInput%specType(ispec)
1898 
1899  y2 = irr*pcvspec(ispec,c2) + pgcspec(xcoord,ispec,c2)*dx2 &
1900  + pgcspec(ycoord,ispec,c2)*dy2 &
1901  + pgcspec(zcoord,ispec,c2)*dz2
1902 
1903  mmr = mmr + y2/pspectype%pMaterial%molw
1904  cpr = cpr + y2*pspectype%pMaterial%spht
1905  END DO ! iSpec
1906 
1907  mmr = 1.0_rfreal/mmr
1908  gcr = mixtperf_r_m(mmr)
1909  gr = mixtperf_g_cpr(cpr,gcr)
1910 #endif
1911 
1912  ur = pcv(cv_mixt_xmom,c2)*irr
1913  vr = pcv(cv_mixt_ymom,c2)*irr
1914  wr = pcv(cv_mixt_zmom,c2)*irr
1915  pr = pdv(dv_mixt_pres,c2)
1916 
1917  rr = rr + pgc(xcoord,grc_mixt_dens,c2)*dx2 &
1918  + pgc(ycoord,grc_mixt_dens,c2)*dy2 &
1919  + pgc(zcoord,grc_mixt_dens,c2)*dz2
1920  ur = ur + pgc(xcoord,grc_mixt_xvel,c2)*dx2 &
1921  + pgc(ycoord,grc_mixt_xvel,c2)*dy2 &
1922  + pgc(zcoord,grc_mixt_xvel,c2)*dz2
1923  vr = vr + pgc(xcoord,grc_mixt_yvel,c2)*dx2 &
1924  + pgc(ycoord,grc_mixt_yvel,c2)*dy2 &
1925  + pgc(zcoord,grc_mixt_yvel,c2)*dz2
1926  wr = wr + pgc(xcoord,grc_mixt_zvel,c2)*dx2 &
1927  + pgc(ycoord,grc_mixt_zvel,c2)*dy2 &
1928  + pgc(zcoord,grc_mixt_zvel,c2)*dz2
1929  pr = pr + pgc(xcoord,grc_mixt_pres,c2)*dx2 &
1930  + pgc(ycoord,grc_mixt_pres,c2)*dy2 &
1931  + pgc(zcoord,grc_mixt_pres,c2)*dz2
1932 
1933  er = mixtperf_eo_dgpuvw(rr,gr,pr,ur,vr,wr)
1934  ar = mixtperf_c_dgp(rr,gr,pr)
1935 
1936  vmr = ur*ur + vr*vr + wr*wr
1937  qr = ur*nx + vr*ny + wr*nz - fs
1938  hr = er + irr*pr
1939 
1940 ! ==============================================================================
1941 ! Compute weights
1942 ! ==============================================================================
1943 
1944  wt1 = sqrt(rr/rl)
1945  wt2 = 1.0_rfreal/(1.0_rfreal + wt1)
1946 
1947 ! ==============================================================================
1948 ! Compute averages. NOTE this average differs from that of the TCP case in
1949 ! that need to recompute gas properties to compute speed of sound. Choose to
1950 ! apply Roe-average also to mass fractions, although there is no analysis to
1951 ! back this up. It is not clear what is the best way to do this.
1952 ! ==============================================================================
1953 
1954  us = wt2*(ul + wt1*ur)
1955  vs = wt2*(vl + wt1*vr)
1956  ws = wt2*(wl + wt1*wr)
1957  qs = wt2*(ql + wt1*qr)
1958  hs = wt2*(hl + wt1*hr)
1959 
1960 #ifdef SPEC
1961  mms = 0.0_rfreal
1962  cps = 0.0_rfreal
1963 
1964  DO ispec = 1,pregion%specInput%nSpecies
1965  pspectype => pregion%specInput%specType(ispec)
1966 
1967  y1 = irl*pcvspec(ispec,c1) + pgcspec(xcoord,ispec,c1)*dx1 &
1968  + pgcspec(ycoord,ispec,c1)*dy1 &
1969  + pgcspec(zcoord,ispec,c1)*dz1
1970 
1971  y2 = irr*pcvspec(ispec,c2) + pgcspec(xcoord,ispec,c2)*dx2 &
1972  + pgcspec(ycoord,ispec,c2)*dy2 &
1973  + pgcspec(zcoord,ispec,c2)*dz2
1974 
1975  ys = wt2*(y1 + wt1*y2)
1976 
1977  mms = mms + ys/pspectype%pMaterial%molw
1978  cps = cps + ys*pspectype%pMaterial%spht
1979  END DO ! iSpec
1980 
1981  mms = 1.0_rfreal/mms
1982  gcs = mixtperf_r_m(mms)
1983  gs = mixtperf_g_cpr(cps,gcs)
1984 #endif
1985 
1986  vms = us*us + vs*vs + ws*ws
1987  as = mixtperf_c_ghovm2(gs,hs,vms)
1988  sl = min(ql-al,qs-as)
1989  sr = max(qr+ar,qs+as)
1990 
1991 ! ==============================================================================
1992 ! Compute fluxes
1993 ! ==============================================================================
1994 
1995  IF ( sl > 0.0_rfreal ) THEN
1996  flx(1) = ql* rl *nm
1997  flx(2) = (ql* rl*ul + pl*nx )*nm
1998  flx(3) = (ql* rl*vl + pl*ny )*nm
1999  flx(4) = (ql* rl*wl + pl*nz )*nm
2000  flx(5) = (ql*(rl*el + pl) + pl*fs)*nm
2001  ELSE IF ( sr < 0.0_rfreal ) THEN
2002  flx(1) = qr* rr *nm
2003  flx(2) = (qr* rr*ur + pr*nx )*nm
2004  flx(3) = (qr* rr*vr + pr*ny )*nm
2005  flx(4) = (qr* rr*wr + pr*nz )*nm
2006  flx(5) = (qr*(rr*er + pr) + pr*fs)*nm
2007  ELSE
2008  sm = (rr*qr*(sr-qr) - rl*ql*(sl-ql) - pr + pl)/(rr*(sr-qr) - rl*(sl-ql))
2009  ps = pl + rl*(ql-sl)*(ql-sm)
2010 
2011  IF ( sm > 0.0_rfreal ) THEN
2012  term = 1.0_rfreal/(sl-sm)
2013 
2014  rs = term* (sl-ql)*rl
2015  rus = term*((sl-ql)*rl*ul + (ps - pl )*nx)
2016  rvs = term*((sl-ql)*rl*vl + (ps - pl )*ny)
2017  rws = term*((sl-ql)*rl*wl + (ps - pl )*nz)
2018  res = term*((sl-ql)*rl*el + (ps*sm - pl*ql) )
2019  ELSE
2020  term = 1.0_rfreal/(sr-sm)
2021 
2022  rs = term* (sr-qr)*rr
2023  rus = term*((sr-qr)*rr*ur + (ps - pr )*nx)
2024  rvs = term*((sr-qr)*rr*vr + (ps - pr )*ny)
2025  rws = term*((sr-qr)*rr*wr + (ps - pr )*nz)
2026  res = term*((sr-qr)*rr*er + (ps*sm - pr*qr) )
2027  END IF ! sm
2028 
2029  flx(1) = sm* rs *nm
2030  flx(2) = (sm* rus + ps*nx )*nm
2031  flx(3) = (sm* rvs + ps*ny )*nm
2032  flx(4) = (sm* rws + ps*nz )*nm
2033  flx(5) = (sm*(res + ps) - ps*fs)*nm
2034  END IF ! sl
2035 
2036 ! ==============================================================================
2037 ! Store mass flux
2038 ! ==============================================================================
2039 
2040  pmf(indmf*ifg) = flx(1)
2041 
2042 ! ==============================================================================
2043 ! Accumulate into residual
2044 ! ==============================================================================
2045 
2046  prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + flx(1)
2047  prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + flx(2)
2048  prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + flx(3)
2049  prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + flx(4)
2050  prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + flx(5)
2051 
2052  prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - flx(1)
2053  prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - flx(2)
2054  prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - flx(3)
2055  prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - flx(4)
2056  prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - flx(5)
2057 
2058 ! ==============================================================================
2059 ! Accumulate into substantial derivative
2060 ! ==============================================================================
2061 
2062  psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + us*pmf(indmf*ifg)
2063  psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vs*pmf(indmf*ifg)
2064  psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + ws*pmf(indmf*ifg)
2065 
2066  psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - us*pmf(indmf*ifg)
2067  psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vs*pmf(indmf*ifg)
2068  psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - ws*pmf(indmf*ifg)
2069  END DO ! ifg
2070 
2071 ! ******************************************************************************
2072 ! Compute fluxes through boundary faces
2073 ! ******************************************************************************
2074 
2075  DO ipatch = 1,pregion%grid%nPatches
2076  ppatch => pregion%patches(ipatch)
2077 
2078  SELECT CASE ( ppatch%spaceOrder )
2079  CASE ( 1 )
2080  CALL rflu_centralfirstpatch(pregion,ppatch)
2081  CASE ( 2 )
2082  CALL rflu_centralsecondpatch(pregion,ppatch)
2083  CASE default
2084  CALL errorstop(global,err_reached_default,__line__)
2085  END SELECT ! pPatch%spaceOrder
2086  END DO ! iPatch
2087 
2088 ! ******************************************************************************
2089 ! End
2090 ! ******************************************************************************
2091 
2092 #ifdef ROCPROF
2093  CALL fprofiler_ends("RFLU::HLLC_ComputeFlux2_MTCP")
2094 #endif
2095 
2096  CALL deregisterfunction(global)
2097 
2098 END SUBROUTINE rflu_hllc_computeflux2_mtcp
2099 
2100 
2101 
2102 
2103 
2104 
2105 
2106 
2107 ! ******************************************************************************
2108 !
2109 ! Purpose: Compute convective fluxes using second-order accurate version of
2110 ! HLLC scheme for thermally and calorically perfect gas.
2111 !
2112 ! Description: None.
2113 !
2114 ! Input:
2115 ! pRegion Pointer to region
2116 !
2117 ! Output: None.
2118 !
2119 ! Notes:
2120 ! 1. Only applicable to thermally and calorically perfect gas because would
2121 ! need to recompute mixture properties at faces depending on face states.
2122 !
2123 ! ******************************************************************************
2124 
2125 SUBROUTINE rflu_hllc_computeflux2_tcp(pRegion)
2126 
2127  USE modinterfaces, ONLY: mixtperf_c_dgp, &
2128  mixtperf_c_ghovm2, &
2129  mixtperf_eo_dgpuvw, &
2130  rflu_centralfirstpatch, &
2131  rflu_centralsecondpatch
2132 
2133  IMPLICIT NONE
2134 
2135 ! ******************************************************************************
2136 ! Definitions and declarations
2137 ! ******************************************************************************
2138 
2139 ! ==============================================================================
2140 ! Arguments
2141 ! ==============================================================================
2142 
2143  TYPE(t_region), POINTER :: pregion
2144 
2145 ! ==============================================================================
2146 ! Locals
2147 ! ==============================================================================
2148 
2149  INTEGER :: c1,c2,ifg,indgs,indmf,indsd,ipatch
2150  REAL(RFREAL) :: al,ar,as,cp,dx,dy,dz,el,er,fs,g,hl,hr,hs,irl,irr,nm, &
2151  nx,ny,nz,ql,qr,qs,pl,pr,ps,rl,rr,res,rs,rus,rvs,rws, &
2152  sl,sm,sr,term,ul,ur,us,vl,vml,vmr,vms,vr,vs,wl,wr,ws, &
2153  wt1,wt2,xc,yc,zc
2154  REAL(RFREAL) :: flx(5)
2155  REAL(RFREAL), DIMENSION(:), POINTER :: pmf
2156  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,prhs,psd
2157  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: pgc
2158  TYPE(t_global), POINTER :: global
2159  TYPE(t_grid), POINTER :: pgrid
2160  TYPE(t_patch), POINTER :: ppatch
2161 
2162 ! ******************************************************************************
2163 ! Start
2164 ! ******************************************************************************
2165 
2166  global => pregion%global
2167 
2168  CALL registerfunction(global,'RFLU_HLLC_ComputeFlux2_TCP',&
2169  'RFLU_ModHLLCFlux.F90')
2170 
2171 #ifdef ROCPROF
2172  CALL fprofiler_begins("RFLU::HLLC_ComputeFlux2_TCP")
2173 #endif
2174 
2175 ! ******************************************************************************
2176 ! Set dimensions and pointers
2177 ! ******************************************************************************
2178 
2179  pgrid => pregion%grid
2180 
2181  indgs = pgrid%indGs
2182 
2183  indmf = pregion%mixtInput%indMfMixt
2184  indsd = pregion%mixtInput%indSd
2185 
2186  pcv => pregion%mixt%cv
2187  pdv => pregion%mixt%dv
2188 
2189  pgc => pregion%mixt%gradCell
2190  pmf => pregion%mixt%mfMixt
2191  prhs => pregion%mixt%rhs
2192  psd => pregion%mixt%sd
2193 
2194  IF ( pregion%mixtInput%gasModel /= gas_model_tcperf ) THEN
2195  CALL errorstop(global,err_gasmodel_invalid,__line__)
2196  END IF ! pRegion%mixtInput%gasModel
2197 
2198  cp = global%refCp
2199  g = global%refGamma
2200 
2201 ! ******************************************************************************
2202 ! Compute fluxes through interior faces
2203 ! ******************************************************************************
2204 
2205  DO ifg = 1,pgrid%nFaces
2206  c1 = pgrid%f2c(1,ifg)
2207  c2 = pgrid%f2c(2,ifg)
2208 
2209 ! ==============================================================================
2210 ! Get face geometry and grid speed
2211 ! ==============================================================================
2212 
2213  nx = pgrid%fn(xcoord,ifg)
2214  ny = pgrid%fn(ycoord,ifg)
2215  nz = pgrid%fn(zcoord,ifg)
2216  nm = pgrid%fn(xyzmag,ifg)
2217 
2218  xc = pgrid%fc(xcoord,ifg)
2219  yc = pgrid%fc(ycoord,ifg)
2220  zc = pgrid%fc(zcoord,ifg)
2221 
2222  fs = pgrid%gs(indgs*ifg)
2223 
2224 ! ==============================================================================
2225 ! Compute left and right states
2226 ! ==============================================================================
2227 
2228 ! ------------------------------------------------------------------------------
2229 ! Left state
2230 ! ------------------------------------------------------------------------------
2231 
2232  rl = pcv(cv_mixt_dens,c1)
2233  irl = 1.0_rfreal/rl
2234 
2235  ul = pcv(cv_mixt_xmom,c1)*irl
2236  vl = pcv(cv_mixt_ymom,c1)*irl
2237  wl = pcv(cv_mixt_zmom,c1)*irl
2238  pl = pdv(dv_mixt_pres,c1)
2239 
2240  dx = xc - pgrid%cofg(xcoord,c1)
2241  dy = yc - pgrid%cofg(ycoord,c1)
2242  dz = zc - pgrid%cofg(zcoord,c1)
2243 
2244  rl = rl + pgc(xcoord,grc_mixt_dens,c1)*dx &
2245  + pgc(ycoord,grc_mixt_dens,c1)*dy &
2246  + pgc(zcoord,grc_mixt_dens,c1)*dz
2247  ul = ul + pgc(xcoord,grc_mixt_xvel,c1)*dx &
2248  + pgc(ycoord,grc_mixt_xvel,c1)*dy &
2249  + pgc(zcoord,grc_mixt_xvel,c1)*dz
2250  vl = vl + pgc(xcoord,grc_mixt_yvel,c1)*dx &
2251  + pgc(ycoord,grc_mixt_yvel,c1)*dy &
2252  + pgc(zcoord,grc_mixt_yvel,c1)*dz
2253  wl = wl + pgc(xcoord,grc_mixt_zvel,c1)*dx &
2254  + pgc(ycoord,grc_mixt_zvel,c1)*dy &
2255  + pgc(zcoord,grc_mixt_zvel,c1)*dz
2256  pl = pl + pgc(xcoord,grc_mixt_pres,c1)*dx &
2257  + pgc(ycoord,grc_mixt_pres,c1)*dy &
2258  + pgc(zcoord,grc_mixt_pres,c1)*dz
2259 
2260  el = mixtperf_eo_dgpuvw(rl,g,pl,ul,vl,wl)
2261  al = mixtperf_c_dgp(rl,g,pl)
2262 
2263  vml = ul*ul + vl*vl + wl*wl
2264  ql = ul*nx + vl*ny + wl*nz - fs
2265  hl = el + irl*pl
2266 
2267 ! ------------------------------------------------------------------------------
2268 ! Right state
2269 ! ------------------------------------------------------------------------------
2270 
2271  rr = pcv(cv_mixt_dens,c2)
2272  irr = 1.0_rfreal/rr
2273 
2274  ur = pcv(cv_mixt_xmom,c2)*irr
2275  vr = pcv(cv_mixt_ymom,c2)*irr
2276  wr = pcv(cv_mixt_zmom,c2)*irr
2277  pr = pdv(dv_mixt_pres,c2)
2278 
2279  dx = xc - pgrid%cofg(xcoord,c2)
2280  dy = yc - pgrid%cofg(ycoord,c2)
2281  dz = zc - pgrid%cofg(zcoord,c2)
2282 
2283  rr = rr + pgc(xcoord,grc_mixt_dens,c2)*dx &
2284  + pgc(ycoord,grc_mixt_dens,c2)*dy &
2285  + pgc(zcoord,grc_mixt_dens,c2)*dz
2286  ur = ur + pgc(xcoord,grc_mixt_xvel,c2)*dx &
2287  + pgc(ycoord,grc_mixt_xvel,c2)*dy &
2288  + pgc(zcoord,grc_mixt_xvel,c2)*dz
2289  vr = vr + pgc(xcoord,grc_mixt_yvel,c2)*dx &
2290  + pgc(ycoord,grc_mixt_yvel,c2)*dy &
2291  + pgc(zcoord,grc_mixt_yvel,c2)*dz
2292  wr = wr + pgc(xcoord,grc_mixt_zvel,c2)*dx &
2293  + pgc(ycoord,grc_mixt_zvel,c2)*dy &
2294  + pgc(zcoord,grc_mixt_zvel,c2)*dz
2295  pr = pr + pgc(xcoord,grc_mixt_pres,c2)*dx &
2296  + pgc(ycoord,grc_mixt_pres,c2)*dy &
2297  + pgc(zcoord,grc_mixt_pres,c2)*dz
2298 
2299  er = mixtperf_eo_dgpuvw(rr,g,pr,ur,vr,wr)
2300  ar = mixtperf_c_dgp(rr,g,pr)
2301 
2302  vmr = ur*ur + vr*vr + wr*wr
2303  qr = ur*nx + vr*ny + wr*nz - fs
2304  hr = er + irr*pr
2305 
2306 ! ==============================================================================
2307 ! Compute weights
2308 ! ==============================================================================
2309 
2310  wt1 = sqrt(rr/rl)
2311  wt2 = 1.0_rfreal/(1.0_rfreal + wt1)
2312 
2313 ! ==============================================================================
2314 ! Compute averages
2315 ! ==============================================================================
2316 
2317  us = wt2*(ul + wt1*ur)
2318  vs = wt2*(vl + wt1*vr)
2319  ws = wt2*(wl + wt1*wr)
2320  qs = wt2*(ql + wt1*qr)
2321  hs = wt2*(hl + wt1*hr)
2322 
2323  vms = us*us + vs*vs + ws*ws
2324  as = mixtperf_c_ghovm2(g,hs,vms)
2325  sl = min(ql-al,qs-as)
2326  sr = max(qr+ar,qs+as)
2327 
2328 ! ==============================================================================
2329 ! Compute fluxes
2330 ! ==============================================================================
2331 
2332  IF ( sl > 0.0_rfreal ) THEN
2333  flx(1) = ql* rl *nm
2334  flx(2) = (ql* rl*ul + pl*nx )*nm
2335  flx(3) = (ql* rl*vl + pl*ny )*nm
2336  flx(4) = (ql* rl*wl + pl*nz )*nm
2337  flx(5) = (ql*(rl*el + pl) + pl*fs)*nm
2338  ELSE IF ( sr < 0.0_rfreal ) THEN
2339  flx(1) = qr* rr *nm
2340  flx(2) = (qr* rr*ur + pr*nx )*nm
2341  flx(3) = (qr* rr*vr + pr*ny )*nm
2342  flx(4) = (qr* rr*wr + pr*nz )*nm
2343  flx(5) = (qr*(rr*er + pr) + pr*fs)*nm
2344  ELSE
2345  sm = (rr*qr*(sr-qr) - rl*ql*(sl-ql) - pr + pl)/(rr*(sr-qr) - rl*(sl-ql))
2346  ps = pl + rl*(ql-sl)*(ql-sm)
2347 
2348  IF ( sm > 0.0_rfreal ) THEN
2349  term = 1.0_rfreal/(sl-sm)
2350 
2351  rs = term* (sl-ql)*rl
2352  rus = term*((sl-ql)*rl*ul + (ps - pl )*nx)
2353  rvs = term*((sl-ql)*rl*vl + (ps - pl )*ny)
2354  rws = term*((sl-ql)*rl*wl + (ps - pl )*nz)
2355  res = term*((sl-ql)*rl*el + (ps*sm - pl*ql) )
2356  ELSE
2357  term = 1.0_rfreal/(sr-sm)
2358 
2359  rs = term* (sr-qr)*rr
2360  rus = term*((sr-qr)*rr*ur + (ps - pr )*nx)
2361  rvs = term*((sr-qr)*rr*vr + (ps - pr )*ny)
2362  rws = term*((sr-qr)*rr*wr + (ps - pr )*nz)
2363  res = term*((sr-qr)*rr*er + (ps*sm - pr*qr) )
2364  END IF ! sm
2365 
2366  flx(1) = sm* rs *nm
2367  flx(2) = (sm* rus + ps*nx )*nm
2368  flx(3) = (sm* rvs + ps*ny )*nm
2369  flx(4) = (sm* rws + ps*nz )*nm
2370  flx(5) = (sm*(res + ps) - ps*fs)*nm
2371  END IF ! sl
2372 
2373 ! ==============================================================================
2374 ! Store mass flux
2375 ! ==============================================================================
2376 
2377  pmf(indmf*ifg) = flx(1)
2378 
2379 ! ==============================================================================
2380 ! Accumulate into residual
2381 ! ==============================================================================
2382 
2383  prhs(cv_mixt_dens,c1) = prhs(cv_mixt_dens,c1) + flx(1)
2384  prhs(cv_mixt_xmom,c1) = prhs(cv_mixt_xmom,c1) + flx(2)
2385  prhs(cv_mixt_ymom,c1) = prhs(cv_mixt_ymom,c1) + flx(3)
2386  prhs(cv_mixt_zmom,c1) = prhs(cv_mixt_zmom,c1) + flx(4)
2387  prhs(cv_mixt_ener,c1) = prhs(cv_mixt_ener,c1) + flx(5)
2388 
2389  prhs(cv_mixt_dens,c2) = prhs(cv_mixt_dens,c2) - flx(1)
2390  prhs(cv_mixt_xmom,c2) = prhs(cv_mixt_xmom,c2) - flx(2)
2391  prhs(cv_mixt_ymom,c2) = prhs(cv_mixt_ymom,c2) - flx(3)
2392  prhs(cv_mixt_zmom,c2) = prhs(cv_mixt_zmom,c2) - flx(4)
2393  prhs(cv_mixt_ener,c2) = prhs(cv_mixt_ener,c2) - flx(5)
2394 
2395 ! ==============================================================================
2396 ! Accumulate into substantial derivative
2397 ! ==============================================================================
2398 
2399  psd(sd_xmom,c1*indsd) = psd(sd_xmom,c1*indsd) + us*pmf(indmf*ifg)
2400  psd(sd_ymom,c1*indsd) = psd(sd_ymom,c1*indsd) + vs*pmf(indmf*ifg)
2401  psd(sd_zmom,c1*indsd) = psd(sd_zmom,c1*indsd) + ws*pmf(indmf*ifg)
2402 
2403  psd(sd_xmom,c2*indsd) = psd(sd_xmom,c2*indsd) - us*pmf(indmf*ifg)
2404  psd(sd_ymom,c2*indsd) = psd(sd_ymom,c2*indsd) - vs*pmf(indmf*ifg)
2405  psd(sd_zmom,c2*indsd) = psd(sd_zmom,c2*indsd) - ws*pmf(indmf*ifg)
2406  END DO ! ifg
2407 
2408 ! ******************************************************************************
2409 ! Compute fluxes through boundary faces
2410 ! ******************************************************************************
2411 
2412  DO ipatch = 1,pregion%grid%nPatches
2413  ppatch => pregion%patches(ipatch)
2414 
2415  SELECT CASE ( ppatch%spaceOrder )
2416  CASE ( 1 )
2417  CALL rflu_centralfirstpatch(pregion,ppatch)
2418  CASE ( 2 )
2419  CALL rflu_centralsecondpatch(pregion,ppatch)
2420  CASE default
2421  CALL errorstop(global,err_reached_default,__line__)
2422  END SELECT ! pPatch%spaceOrder
2423  END DO ! iPatch
2424 
2425 ! ******************************************************************************
2426 ! End
2427 ! ******************************************************************************
2428 
2429 #ifdef ROCPROF
2430  CALL fprofiler_ends("RFLU::HLLC_ComputeFlux2_TCP")
2431 #endif
2432 
2433  CALL deregisterfunction(global)
2434 
2435 END SUBROUTINE rflu_hllc_computeflux2_tcp
2436 
2437 
2438 
2439 
2440 
2441 
2442 
2443 
2444 ! ******************************************************************************
2445 ! End
2446 ! ******************************************************************************
2447 
2448 END MODULE rflu_modhllcflux
2449 
2450 ! ******************************************************************************
2451 !
2452 ! RCS Revision history:
2453 !
2454 ! $Log: RFLU_ModHLLCFlux.F90,v $
2455 ! Revision 1.10 2008/12/06 08:44:22 mtcampbe
2456 ! Updated license.
2457 !
2458 ! Revision 1.9 2008/11/19 22:17:33 mtcampbe
2459 ! Added Illinois Open Source License/Copyright
2460 !
2461 ! Revision 1.8 2006/05/01 22:15:29 haselbac
2462 ! Bug fix: Moved gs out of SPEC declaration, removed debug statements
2463 !
2464 ! Revision 1.7 2006/05/01 21:00:47 haselbac
2465 ! Rewrite for consistency and cleanliness
2466 !
2467 ! Revision 1.6 2006/04/15 16:58:42 haselbac
2468 ! Added capability of running 1st order boundary fluxes with 2nd order volume
2469 ! fluxes
2470 !
2471 ! Revision 1.5 2006/04/07 15:19:19 haselbac
2472 ! Removed tabs
2473 !
2474 ! Revision 1.4 2006/03/30 20:49:40 haselbac
2475 ! Bug fixes for perf gas routines, cosmetics elsewhere
2476 !
2477 ! Revision 1.3 2006/03/26 20:22:02 haselbac
2478 ! Added support for GL model
2479 !
2480 ! Revision 1.2 2005/07/07 22:45:01 haselbac
2481 ! Added profiling calls
2482 !
2483 ! Revision 1.1 2005/05/16 20:36:29 haselbac
2484 ! Initial revision
2485 ! ******************************************************************************
2486 
2487 
2488 
2489 
2490 
2491 
2492 
2493 
2494 
2495 
2496 
2497 
2498 
2499 
2500 
rflu_centralsecondpatch
subroutine rflu_centralsecondpatch(pRegion, pPatch)
Definition: RFLU_CentralSecondPatch.F90:46
mixtperf_r_m
real(rfreal) function mixtperf_r_m(M)
Definition: MixtPerf_R.F90:54
rs
subroutine rs(nm, n, a, w, matz, z, fv1, fv2, ierr)
Definition: arruda_boyce.f90:259
dx
NT dx
Definition: rational_rotation.h:81
Kdtree_d::ExtPoint::term
void term()
Definition: kdtree_d.h:237
max
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
mixtperf_c_dgp
real(rfreal) function mixtperf_c_dgp(D, G, P)
Definition: MixtPerf_C.F90:56
moderror::registerfunction
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
modgrid::t_grid
Definition: ModGrid.F90:56
mixtperf_d_prt
real(rfreal) function mixtperf_d_prt(P, R, T)
Definition: MixtPerf_D.F90:71
sqrt
double sqrt(double d)
Definition: double.h:73
modspecies::t_spec_type
Definition: ModSpecies.F90:64
rflu_modhllcflux::rflu_hllc_computeflux
subroutine, public rflu_hllc_computeflux(pRegion)
Definition: RFLU_ModHLLCFlux.F90:89
mixtperf_c_ghovm2
real(rfreal) function mixtperf_c_ghovm2(G, Ho, Vm2)
Definition: MixtPerf_C.F90:71
mixtperf_r_cpg
real(rfreal) function mixtperf_r_cpg(Cp, G)
Definition: MixtPerf_R.F90:39
mixtgasliq_p
real(rfreal) function mixtgasliq_p(DYl, DYv, DYg, Cl2, Cv2, Cg2, D, Dz, Po, To, Bp, Bt, T)
Definition: MixtGasLiq_P.F90:40
rflu_centralsecondpatch_gl
subroutine rflu_centralsecondpatch_gl(pRegion, pPatch)
Definition: RFLU_CentralSecondPatch_GL.F90:46
rflu_modhllcflux::rflu_hllc_computeflux1_mtcp
subroutine rflu_hllc_computeflux1_mtcp(pRegion)
Definition: RFLU_ModHLLCFlux.F90:561
rflu_centralfirstpatch_gl
subroutine rflu_centralfirstpatch_gl(pRegion, pPatch)
Definition: RFLU_CentralFirstPatch_GL.F90:46
mixtperf_c2_grt
real(rfreal) function mixtperf_c2_grt(G, R, T)
Definition: MixtPerf_C.F90:101
rflu_modhllcflux::rflu_hllc_computeflux1_tcp
subroutine rflu_hllc_computeflux1_tcp(pRegion)
Definition: RFLU_ModHLLCFlux.F90:905
mixtliq_c2_bp
real(rfreal) function mixtliq_c2_bp(Bp)
Definition: MixtLiq_C.F90:54
rflu_modhllcflux::rflu_hllc_computeflux2_tcp
subroutine rflu_hllc_computeflux2_tcp(pRegion)
Definition: RFLU_ModHLLCFlux.F90:2125
rflu_modhllcflux::rflu_hllc_computeflux2_mtcp
subroutine rflu_hllc_computeflux2_mtcp(pRegion)
Definition: RFLU_ModHLLCFlux.F90:1692
Direction_3::dz
RT dz() const
Definition: Direction_3.h:133
min
Vector_n min(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:346
mixtliq_d_dobpppobttto
real(rfreal) function mixtliq_d_dobpppobttto(Dz, Bp, Bt, P, Po, T, To)
Definition: MixtLiq_D.F90:40
mixtgasliq_c
real(rfreal) function mixtgasliq_c(Cvm, D, P, Dl, Dv, Dg, VFl, VFv, VFg, Cl2, Cv2, Cg2, Bl2, Bv2, Bg2)
Definition: MixtGasLiq_C.F90:39
dy
NT dy
Definition: rational_rotation.h:82
mixtperf_eo_dgpuvw
real(rfreal) function mixtperf_eo_dgpuvw(D, G, P, U, V, W)
Definition: MixtPerf_E.F90:40
modglobal::t_global
Definition: ModGlobal.F90:50
rflu_centralfirstpatch
subroutine rflu_centralfirstpatch(pRegion, pPatch)
Definition: RFLU_CentralFirstPatch.F90:46
moderror::errorstop
subroutine errorstop(global, errorCode, errorLine, addMessage)
Definition: ModError.F90:483
rflu_modhllcflux::rflu_hllc_computeflux2_gl
subroutine rflu_hllc_computeflux2_gl(pRegion)
Definition: RFLU_ModHLLCFlux.F90:1181
moderror::deregisterfunction
subroutine deregisterfunction(global)
Definition: ModError.F90:469
rflu_modhllcflux
Definition: RFLU_ModHLLCFlux.F90:39
mixtperf_g_cpr
real(rfreal) function mixtperf_g_cpr(Cp, R)
Definition: MixtPerf_G.F90:39
mixtperf_cv_cpr
real(rfreal) function mixtperf_cv_cpr(Cp, R)
Definition: MixtPerf_Cv.F90:39
rflu_modhllcflux::rflu_hllc_computeflux1_gl
subroutine rflu_hllc_computeflux1_gl(pRegion)
Definition: RFLU_ModHLLCFlux.F90:202
Color::g
unsigned char g() const
Definition: Color.h:69
modbndpatch::t_patch
Definition: ModBndPatch.F90:91