Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PERI_ModHybridDES.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: Suite of rocperi routines associated with hybrid DES.
26 !
27 ! Description: None.
28 !
29 ! Notes: None.
30 !
31 ! ******************************************************************************
32 !
33 ! $Id: PERI_ModHybridDES.F90,v 1.9 2008/12/06 08:44:36 mtcampbe Exp $
34 !
35 ! Copyright: (c) 2004 by the University of Illinois
36 !
37 ! ******************************************************************************
38 
39 MODULE peri_modhybriddes
40 
41  USE moddatatypes
42  USE modglobal, ONLY: t_global
43  USE moddatastruct, ONLY: t_region
44  USE modparameters
45  USE peri_modparameters
46 #ifdef TURB
47  USE turb_modparameters
48 #endif
49  USE moderror
50  USE modmpi
51 
52  IMPLICIT NONE
53 
54  PRIVATE
55 #ifdef RFLO
56  PUBLIC :: peri_rflo_readmean, &
57  peri_comeancorrection
58 #endif
59 
60 ! ******************************************************************************
61 ! Declarations and definitions
62 ! ******************************************************************************
63 
64  CHARACTER(CHRLEN) :: RCSIdentString = &
65  '$RCSfile: PERI_ModHybridDES.F90,v $ $ $'
66 
67 ! ******************************************************************************
68 ! Routines
69 ! ******************************************************************************
70 
71 #ifdef RFLO
72  CONTAINS
73 
74 !******************************************************************************
75 !
76 ! Purpose: read in mean flow solution
77 !
78 ! Description: the following solution formats are supported:
79 ! - RocfloMP ASCII
80 ! - RocfloMP binary.
81 !
82 ! Input: regions = dimensions of all regions.
83 !
84 ! Output: region%levels%peri%cvMean = mean conservative variables (current
85 ! grid level)
86 !
87 ! Notes: solution and grid speeds are read in only for the current grid level;
88 ! solution is also read in for all dummy cells.
89 !
90 !******************************************************************************
91 
92 SUBROUTINE peri_rflo_readmean( regions )
93 
94  USE modinterfaces, ONLY : rflo_readdatafileint, rflo_readdatafilereal
95  IMPLICIT NONE
96 
97 #include "Indexing.h"
98 
99 ! ... parameters
100  TYPE(t_region), POINTER :: regions(:)
101 
102 ! ... loop variables
103  INTEGER :: iReg, i, j, k, n
104 
105 ! ... local variables
106  CHARACTER(2*CHRLEN+17) :: fname
107  CHARACTER(CHRLEN) :: msg
108 
109 #ifdef MPI
110  INTEGER :: status(MPI_STATUS_SIZE)
111 #endif
112  INTEGER :: iLev, iRegFile, ipc, jpc, kpc, nDumCells, iOff, ijOff, ijk
113  INTEGER :: idcbeg, jdcbeg, kdcbeg, idcend, jdcend, kdcend
114  INTEGER :: nDimC, nRvar, intvar(8), errorFlag
115  INTEGER, ALLOCATABLE :: ivar(:,:)
116 
117  REAL(RFREAL) :: rnik
118  REAL(RFREAL), POINTER :: cvMean(:,:)
119  REAL(RFREAL), ALLOCATABLE :: rvar(:,:), cvFile(:,:)
120  LOGICAL :: doread
121 
122  TYPE(t_global), POINTER :: global
123 
124 !******************************************************************************
125 
126  global => regions(1)%global
127 
128  CALL registerfunction( global,'PERI_RFLO_ReadMean',&
129  'PERI_ModHybridDES.F90' )
130 
131 
132 ! check if reading mean flow is applicable ------------------------------------
133 
134  doread = .false.
135 #ifdef TURB
136  DO ireg=1,global%nRegions
137  IF (regions(ireg)%periInput%flowKind /= peri_flow_none) doread = .true.
138  IF ((doread .EQV. .true.) .AND. &
139  (regions(ireg)%mixtInput%turbModel == turb_model_hdessa)) THEN
140  EXIT
141  ELSE
142  doread = .false.
143  ENDIF
144  ENDDO
145 #endif
146  IF (doread .EQV. .false.) goto 999
147 
148 ! allocate fixed-size temporary data arrays -----------------------------------
149 
150  nrvar = 3
151  ALLOCATE( ivar(5,1),stat=errorflag )
152  ALLOCATE( rvar(nrvar,1),stat=errorflag )
153  global%error = errorflag
154  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
155 
156 ! open solution file (only master proc.) --------------------------------------
157 
158  IF (global%myProcid == masterproc) THEN
159 
160  IF (global%solutFormat == format_ascii) THEN
161  WRITE(fname,'(A,1PE11.5)') trim(global%inDir)//trim(global%casename)//'.msola'
162  OPEN(if_solut,file=fname,form='formatted',status='old',iostat=errorflag)
163  ELSE IF (global%solutFormat == format_binary) THEN
164  WRITE(fname,'(A,1PE11.5)') trim(global%inDir)//trim(global%casename)//'.msolb'
165  OPEN(if_solut,file=fname,form='unformatted',status='old',iostat=errorflag)
166  ELSE
167  CALL errorstop( global,err_unknown_format,__line__ )
168  ENDIF
169 
170  global%error = errorflag
171  IF (global%error /= 0) &
172  CALL errorstop( global,err_file_open,__line__,'File: '//trim(fname) )
173 
174  ENDIF ! MASTERPROC
175 
176 ! read & broadcast time and initial residual in file --------------------------
177 
178  IF (global%myProcid == masterproc) THEN
179  CALL rflo_readdatafilereal( global,if_solut,global%solutFormat,nrvar,1,rvar )
180  ENDIF
181 
182 ! read solution data ----------------------------------------------------------
183 
184  DO ireg=1,global%nRegions
185  ilev = regions(ireg)%currLevel
186 
187 ! - read region number and dimensions (only master)
188 
189  IF (global%myProcid == masterproc) THEN
190  CALL rflo_readdatafileint( global,if_solut,global%solutFormat,5,1,ivar )
191  iregfile = ivar(1,1)
192  ipc = ivar(2,1)
193  jpc = ivar(3,1)
194  kpc = ivar(4,1)
195  ndumcells = ivar(5,1)
196 
197 ! --- get dimensions and pointers
198 
199  idcbeg = 1-ndumcells
200  idcend = ipc+ndumcells
201  jdcbeg = 1-ndumcells
202  jdcend = jpc+ndumcells
203  kdcbeg = 1-ndumcells
204  kdcend = kpc+ndumcells
205  ndimc = (idcend-idcbeg+1)*(jdcend-jdcbeg+1)*(kdcend-kdcbeg+1)
206 
207  IF (iregfile /= ireg) &
208  CALL errorstop( global,err_region_number,__line__,'File: '//trim(fname) )
209  IF ((ipc /= regions(ireg)%levels(ilev)%grid%ipc) .OR. &
210  (jpc /= regions(ireg)%levels(ilev)%grid%jpc)) THEN
211  WRITE(msg,1005) ireg,ipc,jpc
212  CALL errorstop( global,err_grid_dimensions,__line__,msg )
213  ENDIF
214  IF (ndumcells /= regions(ireg)%nDumCells) THEN
215  WRITE(msg,1010) ireg,ndumcells,regions(ireg)%nDumCells
216  CALL errorstop( global,err_grid_dumcells,__line__,msg )
217  ENDIF
218  ENDIF
219 
220 ! - master reads & sends data, others receive them
221 
222  IF (global%myProcid == masterproc) THEN
223 
224  ALLOCATE( cvfile(cv_mixt_neqs,ndimc),stat=errorflag )
225  global%error = errorflag
226  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
227 
228  CALL rflo_readdatafilereal( global,if_solut,global%solutFormat, &
229  cv_mixt_neqs,ndimc,cvfile )
230 
231 #ifdef MPI
232  IF (regions(ireg)%procid /= masterproc) THEN
233  intvar(1) = idcbeg
234  intvar(2) = idcend
235  intvar(3) = jdcbeg
236  intvar(4) = jdcend
237  intvar(5) = kdcbeg
238  intvar(6) = kdcend
239  intvar(7) = ndimc
240  intvar(8) = ndumcells
241  CALL mpi_send( intvar,8,mpi_integer,regions(ireg)%procid, &
242  global%nRegions+ireg,global%mpiComm,global%mpierr )
243  IF (global%mpierr /=0 ) CALL errorstop( global,err_mpi_trouble,__line__ )
244 
245  CALL mpi_send( cvfile,cv_mixt_neqs*ndimc,mpi_rfreal, &
246  regions(ireg)%procid,ireg, &
247  global%mpiComm,global%mpierr )
248  IF (global%mpierr /=0 ) CALL errorstop( global,err_mpi_trouble,__line__ )
249  ENDIF
250 #endif
251 
252  ELSE ! not the master
253 
254  IF (regions(ireg)%procid == global%myProcid) THEN
255 #ifdef MPI
256  CALL mpi_recv( intvar,8,mpi_integer,masterproc,global%nRegions+ireg, &
257  global%mpiComm,status,global%mpierr )
258  IF (global%mpierr /=0 ) CALL errorstop( global,err_mpi_trouble,__line__ )
259  idcbeg = intvar(1)
260  idcend = intvar(2)
261  jdcbeg = intvar(3)
262  jdcend = intvar(4)
263  kdcbeg = intvar(5)
264  kdcend = intvar(6)
265  ndimc = intvar(7)
266  ndumcells = intvar(8)
267 #endif
268  ALLOCATE( cvfile(cv_mixt_neqs,ndimc),stat=errorflag )
269  global%error = errorflag
270  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
271 #ifdef MPI
272  CALL mpi_recv( cvfile,cv_mixt_neqs*ndimc,mpi_rfreal,masterproc,ireg, &
273  global%mpiComm,status,global%mpierr )
274  IF (global%mpierr /=0 ) CALL errorstop( global,err_mpi_trouble,__line__ )
275 #endif
276  ENDIF
277 
278  ENDIF
279 
280 ! - copy solution into data structure
281 
282  IF (regions(ireg)%procid == global%myProcid) THEN
283  cvmean => regions(ireg)%levels(ilev)%peri%cvMean
284  cvmean = 0._rfreal
285  ioff = idcend-idcbeg+1
286  ijoff = ioff*(jdcend-jdcbeg+1)
287  n = 0
288  rnik = 1._rfreal/REAL((idcend-idcbeg+1-2*ndumCells)* &
289  (kdcend-kdcbeg+1-2*ndumCells))
290  DO k=kdcbeg+ndumcells,kdcend-ndumcells
291  DO j=jdcbeg,jdcend
292  DO i=idcbeg+ndumcells,idcend-ndumcells
293  n = n + 1
294  ijk = indijk(i,j,k,ioff,ijoff)- &
295  indijk(idcbeg,jdcbeg,kdcbeg,ioff,ijoff) + 1
296  cvmean(j,cv_mixt_dens) = cvmean(j,cv_mixt_dens)+cvfile(1,ijk)
297  cvmean(j,cv_mixt_xmom) = cvmean(j,cv_mixt_xmom)+cvfile(2,ijk)
298  cvmean(j,cv_mixt_ymom) = cvmean(j,cv_mixt_ymom)+cvfile(3,ijk)
299  cvmean(j,cv_mixt_zmom) = cvmean(j,cv_mixt_zmom)+cvfile(4,ijk)
300  cvmean(j,cv_mixt_ener) = cvmean(j,cv_mixt_ener)+cvfile(5,ijk)
301  ENDDO
302  ENDDO
303  ENDDO
304  DO j=jdcbeg,jdcend
305  cvmean(j,cv_mixt_dens) = cvmean(j,cv_mixt_dens)*rnik
306  cvmean(j,cv_mixt_xmom) = cvmean(j,cv_mixt_xmom)*rnik
307  cvmean(j,cv_mixt_ymom) = cvmean(j,cv_mixt_ymom)*rnik
308  cvmean(j,cv_mixt_zmom) = cvmean(j,cv_mixt_zmom)*rnik
309  cvmean(j,cv_mixt_ener) = cvmean(j,cv_mixt_ener)*rnik
310  ENDDO
311  ENDIF ! global%myProcid
312 
313  IF (ALLOCATED(cvfile)) THEN
314  DEALLOCATE( cvfile,stat=errorflag )
315  global%error = errorflag
316  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
317  ENDIF
318  ENDDO ! iReg
319 
320 ! finalize --------------------------------------------------------------------
321 
322  IF (global%myProcid == masterproc) THEN
323  CLOSE(if_solut,iostat=errorflag)
324  global%error = errorflag
325  IF (global%error /= 0) &
326  CALL errorstop( global,err_file_close,__line__,'File: '//trim(fname) )
327  ENDIF
328 
329 999 CONTINUE
330  CALL deregisterfunction( global )
331 
332 1005 FORMAT('Region ',i5,', ipc= ',i6,', jpc= ',i6,', kpc= ',i6,'.')
333 1010 FORMAT('Region ',i5,', # dummy cells=',i2,' but should be= ',i1)
334 
335 END SUBROUTINE peri_rflo_readmean
336 
337 !******************************************************************************
338 !
339 ! Purpose: correct mean flow by RaNS mean flow
340 !
341 ! Description: none.
342 !
343 ! Input: region = data current region
344 !
345 ! Output: region%levels%mixt%cv = conservative variables (current grid level)
346 !
347 ! Notes: tav indexes for rho,ru,rv,rw,pr,rw change if input mixtStatId in
348 ! .input file change.
349 !
350 !******************************************************************************
351 
352 SUBROUTINE peri_comeancorrection( region )
353 
354  USE modinterfaces, ONLY : rflo_getdimensdummy, rflo_getdimensphys, &
355  rflo_getcelloffset, rflo_getnodeoffset, &
356  mixtperf_g_cpr, mixtperf_r_m
357  IMPLICIT NONE
358 
359 #include "Indexing.h"
360 
361 ! ... parameters
362  TYPE(t_region) :: region
363 
364 ! ... loop variables
365  INTEGER :: i, j, k
366 
367 ! ... local variables
368  INTEGER :: iLev, iOff, ijOff, iNOff, ijNOff, ijk, ijkN, ijkNp1
369  INTEGER :: idcbeg, jdcbeg, kdcbeg, idcend, jdcend, kdcend, ibc, iec
370  INTEGER :: ipcbeg, jpcbeg, kpcbeg, ipcend, jpcend, kpcend
371  INTEGER :: indCp, indMol, errorFlag
372  INTEGER, ALLOCATABLE :: jid(:)
373 
374  REAL(RFREAL) :: mm,rgas,cpgas,g, rrho,rho,ru,rv,rw,pr,re, rAvgTim,beta
375  REAL(RFREAL) :: dy1,vistau1,yn1, muel,muet,uv, upl,umi,dudy,ynd,ombeta
376  REAL(RFREAL) :: restau,modtau,vistau,tottau, u,v,w, ampli, alpha,ralpha
377  REAL(RFREAL) :: tavMassflux, cvmMassflux, rnik, fact
378  REAL(RFREAL) :: sndMassFlux, rcvMassFlux, sndAverSize, rcvAverSize
379  REAL(RFREAL), POINTER :: xyz(:,:),cv(:,:),dv(:,:),gv(:,:),tv(:,:)
380  REAL(RFREAL), POINTER :: cvMean(:,:),tav(:,:),ttav(:,:)
381  REAL(RFREAL), ALLOCATABLE :: rhom(:),um(:),vm(:),wm(:),rem(:),dy(:)
382  REAL(RFREAL), ALLOCATABLE :: rhoa(:),ua(:),va(:),wa(:),rea(:)
383  REAL(RFREAL), ALLOCATABLE :: muela(:),mueta(:),uva(:),dstress(:),yc(:),yn(:)
384  LOGICAL :: docorrect
385  TYPE(t_global), POINTER :: global
386 
387 !******************************************************************************
388 
389  global => region%global
390 
391  CALL registerfunction( global,'PERI_CoMeanCorrection',&
392  'PERI_ModHybridDES.F90' )
393 
394 ! check if mean-flow correction is applicable ---------------------------------
395 
396  docorrect = .false.
397 #ifdef TURB
398  IF (region%periInput%flowKind /= peri_flow_none) docorrect = .true.
399  IF ((docorrect .EQV. .true.) .AND. &
400  (region%mixtInput%turbModel == turb_model_hdessa)) THEN
401  ELSE
402  docorrect = .false.
403  ENDIF
404 #endif
405 #ifdef STATS
406  IF ((docorrect .EQV. .true.) .AND. &
407  (global%doStat == active)) THEN
408  ELSE
409  docorrect = .false.
410  ENDIF
411 #else
412  docorrect = .false.
413 #endif
414  IF (docorrect .EQV. .false.) goto 999
415 
416 ! correct mean flow -----------------------------------------------------------
417 
418  ilev = region%currLevel
419 
420  CALL rflo_getdimensdummy( region,ilev,idcbeg,idcend, &
421  jdcbeg,jdcend,kdcbeg,kdcend )
422  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
423  jpcbeg,jpcend,kpcbeg,kpcend )
424  CALL rflo_getcelloffset( region,ilev,ioff,ijoff )
425  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
426 
427  xyz => region%levels(ilev)%grid%xyz
428  cv => region%levels(ilev)%mixt%cv
429  dv => region%levels(ilev)%mixt%dv
430  gv => region%levels(ilev)%mixt%gv
431  tv => region%levels(ilev)%mixt%tv
432  cvmean => region%levels(ilev)%peri%cvMean
433 
434 #ifdef STATS
435 ! compute mean profiles from tav or space-averaged
436 
437  ALLOCATE( rhom(jdcbeg:jdcend),stat=errorflag )
438  global%error = errorflag
439  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
440  ALLOCATE( um(jdcbeg:jdcend),stat=errorflag )
441  global%error = errorflag
442  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
443  ALLOCATE( vm(jdcbeg:jdcend),stat=errorflag )
444  global%error = errorflag
445  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
446  ALLOCATE( wm(jdcbeg:jdcend),stat=errorflag )
447  global%error = errorflag
448  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
449  ALLOCATE( rem(jdcbeg:jdcend),stat=errorflag )
450  global%error = errorflag
451  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
452 
453  ALLOCATE( rhoa(jdcbeg:jdcend),stat=errorflag )
454  global%error = errorflag
455  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
456  ALLOCATE( ua(jdcbeg:jdcend),stat=errorflag )
457  global%error = errorflag
458  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
459  ALLOCATE( va(jdcbeg:jdcend),stat=errorflag )
460  global%error = errorflag
461  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
462  ALLOCATE( wa(jdcbeg:jdcend),stat=errorflag )
463  global%error = errorflag
464  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
465  ALLOCATE( rea(jdcbeg:jdcend),stat=errorflag )
466  global%error = errorflag
467  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
468 
469  ALLOCATE( dy(jdcbeg:jdcend),stat=errorflag )
470  global%error = errorflag
471  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
472  ALLOCATE( yc(jdcbeg:jdcend),stat=errorflag )
473  global%error = errorflag
474  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
475  ALLOCATE( yn(jdcbeg:jdcend+1),stat=errorflag )
476  global%error = errorflag
477  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
478  ALLOCATE( muela(jdcbeg:jdcend),stat=errorflag )
479  global%error = errorflag
480  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
481  ALLOCATE( mueta(jdcbeg:jdcend),stat=errorflag )
482  global%error = errorflag
483  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
484  ALLOCATE( uva(jdcbeg:jdcend),stat=errorflag )
485  global%error = errorflag
486  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
487  ALLOCATE( dstress(jdcbeg:jdcend),stat=errorflag )
488  global%error = errorflag
489  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
490 
491  ALLOCATE( jid(jdcbeg-1:jdcend+1),stat=errorflag )
492  global%error = errorflag
493  IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )
494 
495  tav => region%levels(ilev)%mixt%tav
496  ttav => region%levels(ilev)%turb%tav
497 
498  indcp = region%levels(ilev)%mixt%indCp
499  indmol = region%levels(ilev)%mixt%indMol
500  beta = 1._rfreal
501  ombeta = 1._rfreal-beta
502  alpha = 1._rfreal
503  ralpha = 1._rfreal/alpha
504  rnik = 1._rfreal/REAL((ipcend-ipcbeg+1)*(kpcend-kpcbeg+1))
505 
506  DO j=jdcbeg,jdcend
507 
508  rhom(j) = 0._rfreal
509  um(j) = 0._rfreal
510  vm(j) = 0._rfreal
511  wm(j) = 0._rfreal
512  rem(j) = 0._rfreal
513 
514  DO k=kpcbeg,kpcend
515  DO i=ipcbeg,ipcend
516  ijk = indijk(i,j,k,ioff,ijoff)
517  mm = gv(gv_mixt_mol ,ijk*indmol)
518  rgas = mixtperf_r_m( mm )
519  cpgas = gv(gv_mixt_cp ,ijk*indcp )
520  g = mixtperf_g_cpr( cpgas,rgas )
521 
522  rho = cv(cv_mixt_dens,ijk)
523  ru = cv(cv_mixt_xmom,ijk)
524  rv = cv(cv_mixt_ymom,ijk)
525  rw = cv(cv_mixt_zmom,ijk)
526  re = cv(cv_mixt_ener,ijk)
527 
528  rrho = 1._rfreal/rho
529 
530  rhom(j) = rhom(j) + rho
531  um(j) = um(j) + ru*rrho
532  vm(j) = vm(j) + rv*rrho
533  wm(j) = wm(j) + rw*rrho
534  rem(j) = rem(j) + re
535 
536  ENDDO
537  ENDDO
538  ENDDO
539 
540  IF (global%integrTime < 1.e-13_rfreal) THEN
541  ravgtim = 0._rfreal
542 
543  DO j=jdcbeg,jdcend
544 
545  rhoa(j) = 0._rfreal
546  ua(j) = 0._rfreal
547  va(j) = 0._rfreal
548  wa(j) = 0._rfreal
549  rea(j) = 0._rfreal
550 
551  muela(j) = 0._rfreal
552  mueta(j) = 0._rfreal
553  uva(j) = 0._rfreal
554 
555  DO k=kpcbeg,kpcend
556  DO i=ipcbeg,ipcend
557  ijk = indijk(i,j,k,ioff,ijoff)
558  mm = gv(gv_mixt_mol ,ijk*indmol)
559  rgas = mixtperf_r_m( mm )
560  cpgas = gv(gv_mixt_cp ,ijk*indcp )
561  g = mixtperf_g_cpr( cpgas,rgas )
562 
563 ! rho = tav(1,ijk)*rAvgTim
564 ! ru = rho*tav(2,ijk)*rAvgTim
565 ! rv = rho*tav(3,ijk)*rAvgTim
566 ! rw = 0._RFREAL
567 ! pr = rho*rgas*tav(4,ijk)*rAvgTim
568 ! re = pr/(g-1._RFREAL)+0.5_RFREAL*(ru*ru+rv*rv+rw*rw)/rho
569 
570  rho = cv(cv_mixt_dens,ijk)
571  ru = cv(cv_mixt_xmom,ijk)
572  rv = cv(cv_mixt_ymom,ijk)
573  rw = cv(cv_mixt_zmom,ijk)
574  re = cv(cv_mixt_ener,ijk)
575 
576  rrho = 1._rfreal/rho
577 
578  rhoa(j) = rhoa(j) + rho
579  ua(j) = ua(j) + ru*rrho
580  va(j) = va(j) + rv*rrho
581  wa(j) = wa(j) + rw*rrho
582  rea(j) = rea(j) + re
583 
584 ! ------- for stresses
585 ! muel = global%refVisc
586 ! muet = ttav(1,ijk)*rAvgTim
587 ! uv = tav(9,ijk)*rAvgTim
588 
589  muel = tv(tv_mixt_muel,ijk)
590  muet = tv(tv_mixt_muet,ijk)
591  uv = cv(cv_mixt_xmom,ijk)*cv(cv_mixt_ymom,ijk)*rrho**2
592 
593  muela(j) = muela(j) + muel
594  mueta(j) = mueta(j) + muet
595  uva(j) = uva(j) + uv
596 
597  ENDDO
598  ENDDO
599  ENDDO
600  ELSE ! global%integrTime
601  ravgtim = 1._rfreal/global%integrTime
602 
603  DO j=jdcbeg,jdcend
604 
605  rhoa(j) = 0._rfreal
606  ua(j) = 0._rfreal
607  va(j) = 0._rfreal
608  wa(j) = 0._rfreal
609  rea(j) = 0._rfreal
610 
611  muela(j) = 0._rfreal
612  mueta(j) = 0._rfreal
613  uva(j) = 0._rfreal
614 
615  DO k=kpcbeg,kpcend
616  DO i=ipcbeg,ipcend
617  ijk = indijk(i,j,k,ioff,ijoff)
618  mm = gv(gv_mixt_mol ,ijk*indmol)
619  rgas = mixtperf_r_m( mm )
620  cpgas = gv(gv_mixt_cp ,ijk*indcp )
621  g = mixtperf_g_cpr( cpgas,rgas )
622 
623  rho = tav(1,ijk)*ravgtim
624  ru = rho*tav(2,ijk)*ravgtim
625  rv = rho*tav(3,ijk)*ravgtim
626  rw = 0._rfreal
627  pr = rho*rgas*tav(4,ijk)*ravgtim
628  re = pr/(g-1._rfreal)+0.5_rfreal*(ru*ru+rv*rv+rw*rw)/rho
629 
630 ! rho = cv(CV_MIXT_DENS,ijk)
631 ! ru = cv(CV_MIXT_XMOM,ijk)
632 ! rv = cv(CV_MIXT_YMOM,ijk)
633 ! rw = cv(CV_MIXT_ZMOM,ijk)
634 ! re = cv(CV_MIXT_ENER,ijk)
635 
636  rrho = 1._rfreal/rho
637 
638  rhoa(j) = rhoa(j) + rho
639  ua(j) = ua(j) + ru*rrho
640  va(j) = va(j) + rv*rrho
641  wa(j) = wa(j) + rw*rrho
642  rea(j) = rea(j) + re
643 
644 ! ------- for stresses
645  muel = global%refVisc
646  muet = ttav(1,ijk)*ravgtim
647  uv = tav(9,ijk)*ravgtim
648 
649 ! muel = tv(TV_MIXT_MUEL,ijk)
650 ! muet = tv(TV_MIXT_MUET,ijk)
651 ! uv = cv(CV_MIXT_XMOM,ijk)*cv(CV_MIXT_YMOM,ijk)*rrho**2
652 
653  muela(j) = muela(j) + muel
654  mueta(j) = mueta(j) + muet
655  uva(j) = uva(j) + uv
656 
657  ENDDO
658  ENDDO
659  ENDDO
660  ENDIF ! global%integrTime
661 
662 ! space averaged
663 
664  DO j=jdcbeg,jdcend
665  rhom(j) = rhom(j)*rnik
666  um(j) = um(j)*rnik
667  vm(j) = vm(j)*rnik
668  wm(j) = wm(j)*rnik
669  rem(j) = rem(j)*rnik
670 
671  rhoa(j) = rhoa(j)*rnik
672  ua(j) = ua(j)*rnik
673  va(j) = va(j)*rnik
674  wa(j) = wa(j)*rnik
675  rea(j) = rea(j)*rnik
676 
677  muela(j) = muela(j)*rnik
678  mueta(j) = mueta(j)*rnik
679  uva(j) = uva(j)*rnik
680  ENDDO
681 
682 ! prepare tools for massflux ratio and stress-difference computation
683 
684  tavmassflux = 0._rfreal
685  cvmmassflux = 0._rfreal
686  DO k=1,1
687  DO j=jdcbeg,jdcend
688  DO i=1,1
689  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
690  ijknp1= indijk(i ,j+1,k ,inoff,ijnoff)
691  dy(j) = (xyz(ycoord,ijknp1)-xyz(ycoord,ijkn))
692  yc(j) = 0.5_rfreal*(xyz(ycoord,ijknp1)+xyz(ycoord,ijkn))
693  yn(j) = xyz(ycoord,ijkn)
694  yn(j+1) = xyz(ycoord,ijknp1)
695  ENDDO
696  ENDDO
697  ENDDO
698  IF (region%iRegionGlobal == 1) THEN
699  dy1 = dy(1)
700  vistau1 = muela(1)*2._rfreal*ua(1)/dy1
701  yn1 = yn(1)
702  ENDIF
703 
704 #ifdef MPI
705  CALL mpi_bcast( dy1,1,mpi_rfreal,masterproc,global%mpiComm,global%mpierr )
706  CALL mpi_bcast( vistau1,1,mpi_rfreal,masterproc,global%mpiComm,global%mpierr )
707  CALL mpi_bcast( yn1,1,mpi_rfreal,masterproc,global%mpiComm,global%mpierr )
708 #endif
709 
710  DO j=jpcbeg,jpcend
711  tavmassflux = tavmassflux + dy(j)*rhom(j)*um(j)
712  cvmmassflux = cvmmassflux + dy(j)*cvmean(j,cv_mixt_xmom)
713  ENDDO
714 
715 #ifdef MPI
716  sndmassflux = tavmassflux
717  CALL mpi_allreduce( sndmassflux,rcvmassflux,1,mpi_double_precision,mpi_sum, &
718  global%mpiComm, global%mpierr )
719  tavmassflux = rcvmassflux
720 
721  sndmassflux = cvmmassflux
722  CALL mpi_allreduce( sndmassflux,rcvmassflux,1,mpi_double_precision,mpi_sum, &
723  global%mpiComm, global%mpierr )
724  cvmmassflux = rcvmassflux
725 
726  IF (region%periInput%split(jcoord) == off) THEN
727  sndaversize = rnik
728  CALL mpi_allreduce( sndaversize,rcvaversize,1,mpi_double_precision,mpi_sum, &
729  global%mpiComm, global%mpierr )
730  rnik = rcvaversize
731  ENDIF
732 #endif
733 
734 ! massflux ratio
735 
736  fact = tavmassflux/cvmmassflux
737 
738 ! stress difference
739 
740  DO j=jdcbeg,jdcend
741  jid(j) = j
742  ENDDO
743  jid(jdcbeg-1) = jdcbeg
744  jid(jdcend+1) = jdcend
745 
746  DO j=jdcbeg,jdcend
747 ! dy(j) = MAX( dy(j),dy1 )
748 
749  upl = 0.5_rfreal*(ua(jid(j+1)) + ua(j))
750  umi = 0.5_rfreal*(ua(jid(j-1)) + ua(j))
751 
752  IF (j==jdcbeg .OR. j==jdcend) THEN
753  dudy = 2._rfreal*(upl-umi)/dy(j)
754  ELSE
755  dudy = (upl-umi)/dy(j)
756  ENDIF
757 
758  restau = -rhoa(j)*uva(j)
759  modtau = mueta(j)*dudy
760  vistau = muela(j)*dudy
761 
762  tottau = min( restau+modtau+vistau,vistau1 )
763 
764  IF (tottau < 0._rfreal) tottau = max( restau+modtau+vistau,-vistau1 )
765 
766  ynd = 1._rfreal-(yc(j)-yn1)/global%refLength
767 
768 ! dstress(j) = alpha*( ynd - tottau/vistau1 )/ynd
769 
770  IF (ynd >= 0._rfreal) THEN
771  dstress(j) = alpha*( ynd - tottau/vistau1 )/ &
772  max( restau/vistau1,0.05_rfreal )
773  ELSE
774  dstress(j) = alpha*( tottau/vistau1 - ynd )/ &
775  max(-restau/vistau1,0.05_rfreal )
776  ENDIF
777 
778 ! IF (dstress(j) < 0._RFREAL) dstress(j) = dstress(j)*ralpha
779 
780  IF (abs( ynd ) < 0.1_rfreal) THEN
781  dstress(j) = 0._rfreal
782  ENDIF
783 write(*,100)region%iRegionGlobal,(region%iRegionGlobal-1)*8+j,dstress(j),ynd,tottau,vistau1,tottau/vistau1
784 !write(*,100)region%iRegionGlobal,(region%iRegionGlobal-1)*8+j,restau,modtau,vistau,vistau1,tottau/vistau1
785  ENDDO
786 
787 100 FORMAT( 2i5,5f20.10 )
788 
789 ! correct mean
790 
791  DO k=kdcbeg,kdcend
792  DO j=jdcbeg,jdcend
793  DO i=idcbeg,idcend
794  ijk = indijk(i,j,k,ioff,ijoff)
795 
796  ampli = max( 1._rfreal+dstress(j), 0._rfreal )
797 
798  u = dv(dv_mixt_uvel,ijk)
799  v = dv(dv_mixt_vvel,ijk)
800  w = dv(dv_mixt_wvel,ijk)
801 
802  cv(cv_mixt_dens,ijk) = (cv(cv_mixt_dens,ijk)-rhom(j))*1._rfreal+ &
803  beta*cvmean(j,cv_mixt_dens)+ ombeta*rhom(j)
804 
805  rho = cv(cv_mixt_dens,ijk)
806  rrho = 1._rfreal/cvmean(j,cv_mixt_dens)
807 
808 ! cv(CV_MIXT_XMOM,ijk) = (cv(CV_MIXT_XMOM,ijk)-rum(j))*ampli+ &
809 ! beta*cvMean(j,CV_MIXT_XMOM)+ ombeta*rum(j)
810 
811 ! cv(CV_MIXT_YMOM,ijk) = (cv(CV_MIXT_YMOM,ijk)-rvm(j))*ampli+ &
812 ! beta*cvMean(j,CV_MIXT_YMOM)+ ombeta*rvm(j)
813 
814 ! cv(CV_MIXT_ZMOM,ijk) = (cv(CV_MIXT_ZMOM,ijk)-rwm(j))*ampli+ &
815 ! beta*cvMean(j,CV_MIXT_ZMOM)+ ombeta*rwm(j)
816 
817  cv(cv_mixt_xmom,ijk) = rho*( (u-um(j))*ampli + beta* &
818  cvmean(j,cv_mixt_xmom)*rrho + ombeta*um(j) )
819 
820  cv(cv_mixt_ymom,ijk) = rho*( (v-vm(j))*ampli + beta* &
821  cvmean(j,cv_mixt_ymom)*rrho + ombeta*vm(j) )
822 
823  cv(cv_mixt_zmom,ijk) = rho*( (w-wm(j))*ampli + beta* &
824  cvmean(j,cv_mixt_zmom)*rrho + ombeta*wm(j) )
825 
826  cv(cv_mixt_ener,ijk) = (cv(cv_mixt_ener,ijk)-rem(j))*ampli+ &
827  beta*cvmean(j,cv_mixt_ener)+ ombeta*rem(j)
828 
829 ! mm = gv(GV_MIXT_MOL ,ijk*indMol)
830 ! rgas = MixtPerf_R_M( mm )
831 ! cpgas = gv(GV_MIXT_CP ,ijk*indCp )
832 ! g = MixtPerf_G_CpR( cpgas,rgas )
833 
834 ! rho = cv(CV_MIXT_DENS,ijk)
835 ! ru = cv(CV_MIXT_XMOM,ijk)
836 ! rv = cv(CV_MIXT_YMOM,ijk)
837 ! rw = cv(CV_MIXT_ZMOM,ijk)
838 ! pr = dv(DV_MIXT_PRES,ijk)
839 ! re = pr/(g-1._RFREAL)+0.5_RFREAL*(ru*ru+rv*rv+rw*rw)/rho
840 
841 ! cv(CV_MIXT_ENER,ijk) = (re-rem(j))*1._RFREAL+ &
842 ! beta*cvMean(j,CV_MIXT_ENER)+ ombeta*rem(j)
843 
844  cv(cv_mixt_xmom,ijk) = cv(cv_mixt_xmom,ijk)*fact
845  ENDDO
846  ENDDO
847  ENDDO
848 
849 ! deallocate temporary arrays
850 
851  DEALLOCATE( rhom,stat=errorflag )
852  global%error = errorflag
853  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
854  DEALLOCATE( um,stat=errorflag )
855  global%error = errorflag
856  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
857  DEALLOCATE( vm,stat=errorflag )
858  global%error = errorflag
859  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
860  DEALLOCATE( wm,stat=errorflag )
861  global%error = errorflag
862  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
863  DEALLOCATE( rem,stat=errorflag )
864  global%error = errorflag
865  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
866 
867  DEALLOCATE( rhoa,stat=errorflag )
868  global%error = errorflag
869  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
870  DEALLOCATE( ua,stat=errorflag )
871  global%error = errorflag
872  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
873  DEALLOCATE( va,stat=errorflag )
874  global%error = errorflag
875  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
876  DEALLOCATE( wa,stat=errorflag )
877  global%error = errorflag
878  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
879  DEALLOCATE( rea,stat=errorflag )
880  global%error = errorflag
881  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
882 
883  DEALLOCATE( dy,stat=errorflag )
884  global%error = errorflag
885  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
886  DEALLOCATE( yc,stat=errorflag )
887  global%error = errorflag
888  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
889  DEALLOCATE( yn,stat=errorflag )
890  global%error = errorflag
891  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
892  DEALLOCATE( muela,stat=errorflag )
893  global%error = errorflag
894  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
895  DEALLOCATE( mueta,stat=errorflag )
896  global%error = errorflag
897  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
898  DEALLOCATE( uva,stat=errorflag )
899  global%error = errorflag
900  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
901  DEALLOCATE( dstress,stat=errorflag )
902  global%error = errorflag
903  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
904 
905  DEALLOCATE( jid,stat=errorflag )
906  global%error = errorflag
907  IF (global%error /= 0) CALL errorstop( global,err_deallocate,__line__ )
908 
909 ! compute derived values for the mixture in case of dual time-stepping
910 
911  IF (global%flowType == flow_unsteady) THEN
912 ! IF (global%solverType == SOLV_IMPLICIT) THEN ! dual TST
913  ibc = indijk(idcbeg,jdcbeg,kdcbeg,ioff,ijoff)
914  iec = indijk(idcend,jdcend,kdcend,ioff,ijoff)
915  CALL mixtureproperties( region,ibc,iec,.true. )
916 ! ENDIF
917  ENDIF
918 #endif
919 
920 ! finalize --------------------------------------------------------------------
921 
922 999 CONTINUE
923  CALL deregisterfunction( global )
924 
925 END SUBROUTINE peri_comeancorrection
926 #endif
927 
928 ! ******************************************************************************
929 ! End
930 ! ******************************************************************************
931 
932 END MODULE peri_modhybriddes
933 
934 ! ******************************************************************************
935 !
936 ! RCS Revision history:
937 !
938 ! $Log: PERI_ModHybridDES.F90,v $
939 ! Revision 1.9 2008/12/06 08:44:36 mtcampbe
940 ! Updated license.
941 !
942 ! Revision 1.8 2008/11/19 22:17:49 mtcampbe
943 ! Added Illinois Open Source License/Copyright
944 !
945 ! Revision 1.7 2005/04/21 00:46:19 wasistho
946 ! modified PERI_RFLO_ReadMean
947 !
948 ! Revision 1.6 2005/04/06 02:39:16 wasistho
949 ! modified dstress in PERI_CoMeanCorrection
950 !
951 ! Revision 1.5 2005/04/06 01:53:36 wasistho
952 ! added amplified fluctuations
953 !
954 ! Revision 1.4 2005/03/31 17:03:34 wasistho
955 ! apply space averaged i.o. time averaged and computed dv,gv for dual TST
956 !
957 ! Revision 1.3 2005/03/10 02:02:09 wasistho
958 ! bug fixed for readMean and meanCorrection
959 !
960 ! Revision 1.2 2005/03/07 18:27:56 wasistho
961 ! insert ifdef STATS for tav
962 !
963 ! Revision 1.1 2005/03/07 05:08:13 wasistho
964 ! install hybrid DESSA turbulence model
965 !
966 !
967 ! ******************************************************************************
968 
969 
970 
971 
972 
973 
974 
975 
mixtperf_r_m
real(rfreal) function mixtperf_r_m(M)
Definition: MixtPerf_R.F90:54
max
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
moderror::registerfunction
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
rflo_readdatafileint
subroutine rflo_readdatafileint(global, fileId, form, nDim1, nDim2, ivar)
Definition: RFLO_ReadDataFile.F90:49
turb_modparameters
Definition: TURB_ModParameters.F90:39
rflo_getnodeoffset
subroutine rflo_getnodeoffset(region, iLev, iNodeOffset, ijNodeOffset)
Definition: RFLO_GetDimensions.F90:267
peri_modhybriddes
Definition: PERI_ModHybridDES.F90:39
rflo_getdimensdummy
subroutine rflo_getdimensdummy(region, iLev, idcbeg, idcend, jdcbeg, jdcend, kdcbeg, kdcend)
Definition: RFLO_GetDimensions.F90:67
rflo_getcelloffset
subroutine rflo_getcelloffset(region, iLev, iCellOffset, ijCellOffset)
Definition: RFLO_GetDimensions.F90:235
form
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE form
Definition: RFLO_ReadDataFileInt.h:26
peri_modparameters
Definition: PERI_ModParameters.F90:39
min
Vector_n min(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:346
modglobal::t_global
Definition: ModGlobal.F90:50
moderror::errorstop
subroutine errorstop(global, errorCode, errorLine, addMessage)
Definition: ModError.F90:483
mixtureproperties
subroutine mixtureproperties(region, inBeg, inEnd, gasUpdate)
Definition: MixtureProperties.F90:48
moderror::deregisterfunction
subroutine deregisterfunction(global)
Definition: ModError.F90:469
mixtperf_g_cpr
real(rfreal) function mixtperf_g_cpr(Cp, R)
Definition: MixtPerf_G.F90:39
rflo_getdimensphys
subroutine rflo_getdimensphys(region, iLev, ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend)
Definition: RFLO_GetDimensions.F90:109
rflo_readdatafilereal
subroutine rflo_readdatafilereal(global, fileId, form, nDim1, nDim2, var)
Definition: RFLO_ReadDataFile.F90:98