Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLU_ModOLES.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 routines for optimal LES computations.
26 !
27 ! Description: None.
28 !
29 ! Notes: None.
30 !
31 !*******************************************************************************
32 !
33 ! $Id: RFLU_ModOLES.F90,v 1.10 2008/12/06 08:44:23 mtcampbe Exp $
34 !
35 ! Copyright: (c) 2002 by the University of Illinois
36 !
37 !*******************************************************************************
38 
39 MODULE rflu_modoles
40 
41  USE modglobal, ONLY: t_global
42  USE modparameters
43  USE moddatatypes
44  USE moderror
45  USE modtools, ONLY: makenonzero
46  USE modbndpatch, ONLY: t_patch
47  USE moddatastruct, ONLY: t_region
48  USE modgrid, ONLY: t_grid
49  USE modsortsearch
50  USE modtools
51  USE modmpi
52 
53  USE rflu_modoctree
54 
55  IMPLICIT NONE
56 
57  PRIVATE
58  PUBLIC :: rflu_createstencilsweightsoles, &
59  rflu_createintegralsoles, &
60  rflu_buildstencilsoles, &
61  rflu_findprototypefacesoles, &
62  rflu_computegeometrictermsoles, &
63  rflu_buildsymmetrymapsoles, &
64  rflu_enforcesymmetryoles, &
65  rflu_definecorrelation220, &
66  rflu_definecorrelation221, &
67  rflu_definecorrelation32, &
68  rflu_definecorrelation430, &
69  rflu_definecorrelation431, &
70  rflu_definecorrelation432, &
71  rflu_definecorrelation540, &
72  rflu_definecorrelation541, &
73  rflu_definecorrelation542, &
74  rflu_mapk2ij, &
75  rflu_mapl2ijk, &
76  rflu_mapm2ijkl, &
77  rflu_geti1posoles, &
78  rflu_geti4posoles, &
79  rflu_getlposoles, &
80  rflu_getlposinvoles, &
81  rflu_getqposoles, &
82  rflu_getqposinvoles, &
83  rflu_destroystencilsweightsoles, &
84  rflu_allocatedcuhrearrays, &
85  rflu_deallocatedcuhrearrays, &
86  rflu_setmapfunnz2funcorr22, &
87  rflu_setmapfunnz2funcorr32, &
88  rflu_setmapfunnz2funcorr43, &
89  rflu_setmapfunnz2funcorr44
90 
91  SAVE
92 
93 ! ******************************************************************************
94 ! Declarations and definitions
95 ! ******************************************************************************
96 
97  CHARACTER(CHRLEN) :: &
98  RCSIdentString = '$RCSfile: RFLU_ModOLES.F90,v $ $Revision: 1.10 $'
99  TYPE(t_grid), POINTER :: pGrid
100 
101 ! ==============================================================================
102 ! Variables and parameters used by DCUHRE subroutine
103 ! ==============================================================================
104 
105  INTEGER, PUBLIC :: maxCalls,nDim,nEval,workArraySize,workArraySizeNew
106  INTEGER, PARAMETER, PUBLIC :: DCUHRE_LOOP_LIMIT = 5, &
107  MAX_CALLS_FACTOR = 5, &
108  MAX_CALLS_LIMIT = 50000000, &
109  MAX_CALLS_START = 100000, &
110  MIN_CALLS = 0
111  INTEGER, PUBLIC :: dummy(1)
112  INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: symMapI45OLES
113 
114  REAL(RFREAL), PUBLIC :: errAbsReq,errRelReq
115  REAL(RFREAL), DIMENSION(:), ALLOCATABLE, PUBLIC :: lowLim,errAbsEst, &
116  uppLim,integral, &
117  integralNZ,workArray
118 
119 ! ==============================================================================
120 ! Variables and parameters used in integration routines
121 ! ==============================================================================
122 
123  INTEGER, PUBLIC :: nzLoc,nzSgn
124  INTEGER, DIMENSION(3,3), PARAMETER, PUBLIC :: mapSurf2Vol2 = &
125  RESHAPE((/6,4,4,4,6,5,5,5,6/),(/3,3/))
126  INTEGER, DIMENSION(3,6), PARAMETER, PUBLIC :: mapSurf2Vol3 = &
127  RESHAPE((/9,7,7,7,9,8,8,8,9,9,7,7,7,9,8,8,8,9/),(/3,6/))
128  INTEGER, DIMENSION(3,3), PARAMETER, PUBLIC :: kd = & ! Kronecker delta
129  RESHAPE((/1,0,0,0,1,0,0,0,1/),(/3,3/))
130 
131  INTEGER, PUBLIC :: mapFunNZ2FunCorr22(9),mapFunNZ2FunCorr32(9), &
132  mapFunNZ2FunCorr43(27),mapFunNZ2FunCorr44(81)
133 
134  REAL(RFREAL), PUBLIC :: nzMag,nzVal
135  REAL(RFREAL), PARAMETER, PUBLIC :: CONST_KOLMOGOROV = 2.0_RFREAL
136 
137 
138 ! ******************************************************************************
139 ! Routines
140 ! ******************************************************************************
141 
142  CONTAINS
143 
144 
145 ! ******************************************************************************
146 ! Create geometry
147 ! ******************************************************************************
148 
149  SUBROUTINE rflu_createstencilsweightsoles(pRegion)
150 
151  IMPLICIT NONE
152 
153 ! --- parameters
154 
155  TYPE(t_region), POINTER :: pregion
156 
157 ! --- local variables
158 
159  INTEGER :: errorflag,ncells,nfaces
160  TYPE(t_grid), POINTER :: pgrid
161  TYPE(t_global), POINTER :: global
162 
163 ! ==============================================================================
164 ! Start
165 ! ==============================================================================
166 
167  global => pregion%global
168 
169  CALL registerfunction(global,'RFLU_CreateStencilsWeightsOLES',&
170  'RFLU_ModOLES.F90')
171 
172  IF ( global%myProcid == masterproc .AND. &
173  global%verbLevel > verbose_none ) THEN
174  WRITE(stdout,'(A,1X,A)') solver_name
175  WRITE(stdout,'(A,1X,A)') solver_name,'Creating optimal LES '// &
176  'stencils, weights, and mappings...'
177  END IF ! global%verbLevel
178 
179 ! ==============================================================================
180 ! Set grid pointer
181 ! ==============================================================================
182 
183  pgrid => pregion%grid
184 
185 ! ==============================================================================
186 ! Allocate memory for face stencils, cell limits
187 ! ==============================================================================
188 
189 ! TEMPORARY: hard-code stencil
190  ncells = 2
191 ! END TEMPORARY
192 
193  ALLOCATE(pgrid%fsOLES(ncells,pgrid%nFaces),stat=errorflag)
194  global%error = errorflag
195  IF ( global%error /= err_none ) THEN
196  CALL errorstop(global,err_allocate,__line__,'region%grid%fsOLES')
197  END IF ! global%error
198 
199  pgrid%fsOLES(:,:) = 0
200 
201  ALLOCATE(pgrid%intLimOLES(int_lim_low:int_lim_upp,xcoord:zcoord, &
202  pgrid%nCellsTot),stat=errorflag)
203  global%error = errorflag
204  IF ( global%error /= err_none ) THEN
205  CALL errorstop(global,err_allocate,__line__,'region%grid%intLimOLES')
206  END IF ! global%error
207 
208  pgrid%intLimOLES(int_lim_low,:,:) = huge(1.0_rfreal)
209  pgrid%intLimOLES(int_lim_upp,:,:) = -huge(1.0_rfreal)
210 
211  ALLOCATE(pgrid%rhoOLES(pgrid%nFaces),stat=errorflag)
212  global%error = errorflag
213  IF ( global%error /= err_none ) THEN
214  CALL errorstop(global,err_allocate,__line__,'region%grid%intLimOLES')
215  END IF ! global%error
216 
217  pgrid%rhoOLES(:) = huge(1.0_rfreal)
218 
219  ALLOCATE(pgrid%fp2fOLES(3),stat=errorflag)
220  global%error = errorflag
221  IF ( global%error /= err_none ) THEN
222  CALL errorstop(global,err_allocate,__line__,'region%grid%fp2fOLES')
223  END IF ! global%error
224 
225  ALLOCATE(pgrid%f2fpOLES(pgrid%nFaces),stat=errorflag)
226  global%error = errorflag
227  IF ( global%error /= err_none ) THEN
228  CALL errorstop(global,err_allocate,__line__,'region%grid%f2fpOLES')
229  END IF ! global%error
230 
231 ! ==============================================================================
232 ! Weights, NOTE store only for prototype faces
233 ! ==============================================================================
234 
235  nfaces = 3
236 
237  ALLOCATE(pgrid%wtLinOLES(3,3,ncells,nfaces),stat=errorflag)
238  global%error = errorflag
239  IF ( global%error /= err_none ) THEN
240  CALL errorstop(global,err_allocate,__line__,'region%grid%wtLinOLES')
241  END IF ! global%error
242 
243  ALLOCATE(pgrid%wtQuadOLES(3,3,3,ncells,ncells,nfaces),stat=errorflag)
244  global%error = errorflag
245  IF ( global%error /= err_none ) THEN
246  CALL errorstop(global,err_allocate,__line__,'region%grid%wQuadOLES')
247  END IF ! global%error
248 
249  pgrid%wtLinOLES(:,:,:,:) = REAL(crazy_value_int,kind=rfreal)
250  pgrid%wtQuadOLES(:,:,:,:,:,:) = REAL(crazy_value_int,kind=rfreal)
251 
252 ! ==============================================================================
253 ! Symmetry maps for integrals
254 ! ==============================================================================
255 
256  ALLOCATE(symmapi45oles(9*ncells*ncells),stat=errorflag)
257  global%error = errorflag
258  IF ( global%error /= err_none ) THEN
259  CALL errorstop(global,err_allocate,__line__,'symMapI45OLES')
260  END IF ! global%error
261 
262 ! ==============================================================================
263 ! End
264 ! ==============================================================================
265 
266  IF ( global%myProcid == masterproc .AND. &
267  global%verbLevel > verbose_none ) THEN
268  WRITE(stdout,'(A,1X,A)') solver_name,'Creating optimal LES '// &
269  'stencils, weights, and mappings done.'
270  END IF ! global%verbLevel
271 
272  CALL deregisterfunction(global)
273 
274  END SUBROUTINE rflu_createstencilsweightsoles
275 
276 
277 
278 ! ******************************************************************************
279 ! Create integrals. NOTE initialize to crazy values to be able to check
280 ! whether values are entered correctly.
281 ! ******************************************************************************
282 
283  SUBROUTINE rflu_createintegralsoles(pRegion)
284 
285  IMPLICIT NONE
286 
287 ! --- parameters
288 
289  TYPE(t_region), POINTER :: pregion
290 
291 ! --- local variables
292 
293  INTEGER :: errorflag,ncells,ncols,nfaces,nrows
294  TYPE(t_grid), POINTER :: pgrid
295  TYPE(t_global), POINTER :: global
296 
297 ! ==============================================================================
298 ! Start
299 ! ==============================================================================
300 
301  global => pregion%global
302 
303  CALL registerfunction(global,'RFLU_CreateIntegralsOLES',&
304  'RFLU_ModOLES.F90')
305 
306  IF ( global%myProcid == masterproc .AND. &
307  global%verbLevel > verbose_none ) THEN
308  WRITE(stdout,'(A,1X,A)') solver_name
309  WRITE(stdout,'(A,1X,A)') solver_name,'Creating optimal LES '// &
310  'integrals...'
311  END IF ! global%verbLevel
312 
313 ! ==============================================================================
314 ! Set grid pointer
315 ! ==============================================================================
316 
317  pgrid => pregion%grid
318 
319 ! ==============================================================================
320 ! Allocate memory for integrals, NOTE store only for prototype faces
321 ! ==============================================================================
322 
323  nfaces = 3
324  ncells = SIZE(pgrid%fsOLES,1)
325 
326 ! ------------------------------------------------------------------------------
327 ! Integral 1
328 ! ------------------------------------------------------------------------------
329 
330  ALLOCATE(pgrid%int1OLES(3,nfaces,3*ncells),stat=errorflag)
331  global%error = errorflag
332  IF ( global%error /= err_none ) THEN
333  CALL errorstop(global,err_allocate,__line__,'region%grid%int1OLES')
334  END IF ! global%error
335 
336  pgrid%int1OLES(:,:,:) = REAL(crazy_value_int,kind=rfreal)
337 
338 ! ------------------------------------------------------------------------------
339 ! Integral 2
340 ! ------------------------------------------------------------------------------
341 
342  ALLOCATE(pgrid%int20OLES(nfaces,3*ncells,3*ncells),stat=errorflag)
343  global%error = errorflag
344  IF ( global%error /= err_none ) THEN
345  CALL errorstop(global,err_allocate,__line__,'region%grid%int20OLES')
346  END IF ! global%error
347 
348  ALLOCATE(pgrid%int21OLES(nfaces,3*ncells,3*ncells),stat=errorflag)
349  global%error = errorflag
350  IF ( global%error /= err_none ) THEN
351  CALL errorstop(global,err_allocate,__line__,'region%grid%int21OLES')
352  END IF ! global%error
353 
354  pgrid%int20OLES(:,:,:) = REAL(crazy_value_int,kind=rfreal)
355  pgrid%int21OLES(:,:,:) = REAL(crazy_value_int,kind=rfreal)
356 
357 
358 ! ------------------------------------------------------------------------------
359 ! Integral 3
360 ! ------------------------------------------------------------------------------
361 
362  ALLOCATE(pgrid%int31OLES(nfaces,3*ncells,9*ncells*ncells), &
363  stat=errorflag)
364  global%error = errorflag
365  IF ( global%error /= err_none ) THEN
366  CALL errorstop(global,err_allocate,__line__,'region%grid%int31OLES')
367  END IF ! global%error
368 
369  ALLOCATE(pgrid%int32OLES(nfaces,9*ncells*ncells,3*ncells), &
370  stat=errorflag)
371  global%error = errorflag
372  IF ( global%error /= err_none ) THEN
373  CALL errorstop(global,err_allocate,__line__,'region%grid%int32OLES')
374  END IF ! global%error
375 
376  pgrid%int31OLES(:,:,:) = REAL(crazy_value_int,kind=rfreal)
377  pgrid%int32OLES(:,:,:) = REAL(crazy_value_int,kind=rfreal)
378 
379 ! ------------------------------------------------------------------------------
380 ! Integral 4
381 ! ------------------------------------------------------------------------------
382 
383  ALLOCATE(pgrid%int40OLES(3,nfaces,9*ncells*ncells),stat=errorflag)
384  global%error = errorflag
385  IF ( global%error /= err_none ) THEN
386  CALL errorstop(global,err_allocate,__line__,'region%grid%int40OLES')
387  END IF ! global%error
388 
389  ALLOCATE(pgrid%int41OLES(3,nfaces,9*ncells*ncells),stat=errorflag)
390  global%error = errorflag
391  IF ( global%error /= err_none ) THEN
392  CALL errorstop(global,err_allocate,__line__,'region%grid%int41OLES')
393  END IF ! global%error
394 
395  ALLOCATE(pgrid%int42OLES(3,nfaces,9*ncells*ncells),stat=errorflag)
396  global%error = errorflag
397  IF ( global%error /= err_none ) THEN
398  CALL errorstop(global,err_allocate,__line__,'region%grid%int42OLES')
399  END IF ! global%error
400 
401  pgrid%int40OLES(:,:,:) = REAL(crazy_value_int,kind=rfreal)
402  pgrid%int41OLES(:,:,:) = REAL(crazy_value_int,kind=rfreal)
403  pgrid%int42OLES(:,:,:) = REAL(CRAZY_VALUE_INT,KIND=RFREAL)
404 
405 ! ------------------------------------------------------------------------------
406 ! Integral 5
407 ! ------------------------------------------------------------------------------
408 
409  ALLOCATE(pgrid%int50OLES(nfaces,9*ncells*ncells,9*ncells*ncells), &
410  stat=errorflag)
411  global%error = errorflag
412  IF ( global%error /= err_none ) THEN
413  CALL errorstop(global,err_allocate,__line__,'region%grid%int50OLES')
414  END IF ! global%error
415 
416  ALLOCATE(pgrid%int51OLES(nfaces,9*ncells*ncells,9*ncells*ncells), &
417  stat=errorflag)
418  global%error = errorflag
419  IF ( global%error /= err_none ) THEN
420  CALL errorstop(global,err_allocate,__line__,'region%grid%int51OLES')
421  END IF ! global%error
422 
423  ALLOCATE(pgrid%int52OLES(nfaces,9*ncells*ncells,9*ncells*ncells), &
424  stat=errorflag)
425  global%error = errorflag
426  IF ( global%error /= err_none ) THEN
427  CALL errorstop(global,err_allocate,__line__,'region%grid%int52OLES')
428  END IF ! global%error
429 
430  pgrid%int50OLES(:,:,:) = REAL(crazy_value_int,kind=rfreal)
431  pgrid%int51OLES(:,:,:) = REAL(crazy_value_int,kind=rfreal)
432  pgrid%int52OLES(:,:,:) = REAL(CRAZY_VALUE_INT,KIND=RFREAL)
433 
434 ! ------------------------------------------------------------------------------
435 ! LHS and RHS of linear system
436 ! ------------------------------------------------------------------------------
437 
438  nrows = 3*ncells*(1 + 3*ncells)
439  ncols = nrows
440 
441  ALLOCATE(pgrid%lhsOLES(nrows,ncols),stat=errorflag)
442  global%error = errorflag
443  IF ( global%error /= err_none ) THEN
444  CALL errorstop(global,err_allocate,__line__,'region%grid%lhsOLES')
445  END IF ! global%error
446 
447  ALLOCATE(pgrid%lhsInvOLES(ncols,nrows),stat=errorflag)
448  global%error = errorflag
449  IF ( global%error /= err_none ) THEN
450  CALL errorstop(global,err_allocate,__line__,'region%grid%lhsInvOLES')
451  END IF ! global%error
452 
453  ALLOCATE(pgrid%rhsOLES(nrows),stat=errorflag)
454  global%error = errorflag
455  IF ( global%error /= err_none ) THEN
456  CALL errorstop(global,err_allocate,__line__,'region%grid%rhsOLES')
457  END IF ! global%error
458 
459  pgrid%lhsOLES(:,:) = REAL(CRAZY_VALUE_INT,KIND=RFREAL)
460  pgrid%rhsOLES(:) = REAL(crazy_value_int,kind=rfreal)
461 
462 ! ==============================================================================
463 ! End
464 ! ==============================================================================
465 
466  IF ( global%myProcid == masterproc .AND. &
467  global%verbLevel > verbose_none ) THEN
468  WRITE(stdout,'(A,1X,A)') solver_name,'Creating optimal LES '// &
469  'integrals done.'
470  END IF ! global%verbLevel
471 
472  CALL deregisterfunction(global)
473 
474  END SUBROUTINE rflu_createintegralsoles
475 
476 
477 
478 
479 ! ******************************************************************************
480 ! Build stencil
481 ! ******************************************************************************
482 
483  SUBROUTINE rflu_buildstencilsoles(pRegion)
484 
485  IMPLICIT NONE
486 
487 ! --- parameters
488 
489  TYPE(t_region), POINTER :: pregion
490 
491 ! --- local variables
492 
493  INTEGER :: errorflag,ic,ifc,is,loc,locoffset1,locoffset2,stencilsize, &
494  stencilsizetemp
495  INTEGER, DIMENSION(:), ALLOCATABLE :: stenciltemp
496  REAL(RFREAL) :: delfrac,dotp,vecm,xdel,xmax,xmin,ydel,ymax,ymin, &
497  zdel,zmax,zmin
498  REAL(RFREAL) :: dummy(1),vec(3)
499  REAL(RFREAL), DIMENSION(:), ALLOCATABLE :: dist
500  TYPE(t_grid), POINTER :: pgrid
501  TYPE(t_global), POINTER :: global
502 
503 ! ==============================================================================
504 ! Start
505 ! ==============================================================================
506 
507  global => pregion%global
508 
509  CALL registerfunction(global,'RFLU_BuildStencilsOLES',&
510  'RFLU_ModOLES.F90')
511 
512  IF ( global%myProcid == masterproc .AND. &
513  global%verbLevel > verbose_none ) THEN
514  WRITE(stdout,'(A)') solver_name
515  WRITE(stdout,'(A,1X,A)') solver_name,'Building optimal LES stencils...'
516  END IF ! global%verbLevel
517 
518 ! ==============================================================================
519 ! Set grid pointer and constants
520 ! ==============================================================================
521 
522  pgrid => pregion%grid
523 
524  delfrac = 0.01_rfreal
525 
526 ! ==============================================================================
527 ! Build face stencils
528 ! ==============================================================================
529 
530  stencilsize = SIZE(pgrid%fsOLES,1)
531 
532  IF ( global%myProcid == masterproc .AND. &
533  global%verbLevel > verbose_none ) THEN
534  WRITE(stdout,'(A,3X,A)') solver_name,'Building face stencil support...'
535  WRITE(stdout,'(A,5X,A,1X,I2)') solver_name,'Stencil width:',stencilsize
536  END IF ! global%verbLevel
537 
538 ! ------------------------------------------------------------------------------
539 ! Stencil only contains cells which straddle face, Octree not necessary
540 ! ------------------------------------------------------------------------------
541 
542  IF ( stencilsize == 2 ) THEN ! only need cells which straddle face
543  DO ifc = 1,pgrid%nFaces
544  pgrid%fsOLES(1,ifc) = pgrid%f2c(1,ifc)
545  pgrid%fsOLES(2,ifc) = pgrid%f2c(2,ifc)
546  END DO ! ifc
547 
548 ! ------------------------------------------------------------------------------
549 ! Stencil contains more cells, use Octree
550 ! ------------------------------------------------------------------------------
551 
552  ELSE IF ( stencilsize > 2 ) THEN
553  IF ( stencilsize == 4 ) THEN
554  stencilsizetemp = 20 ! to make sure to capture distance-two cells
555  ELSE IF ( stencilsize == 6 ) THEN
556  stencilsizetemp = 102 ! to make sure to capture distance-three cells
557  ELSE IF ( stencilsize == 8 ) THEN
558  stencilsizetemp = 296 ! to make sure to capture distance-four cells
559  ELSE
560  CALL errorstop(global,err_reached_default,__line__)
561  END IF ! stencilSize
562 
563 ! ----- Allocate memory -------------------------------------------------------
564 
565  ALLOCATE(stenciltemp(stencilsizetemp),stat=errorflag)
566  global%error = errorflag
567  IF ( global%error /= err_none ) THEN
568  CALL errorstop(global,err_allocate,__line__,'stencilTemp')
569  END IF ! global%error
570 
571  ALLOCATE(dist(stencilsizetemp),stat=errorflag)
572  global%error = errorflag
573  IF ( global%error /= err_none ) THEN
574  CALL errorstop(global,err_allocate,__line__,'dist')
575  END IF ! global%error
576 
577 ! ----- Determine bounding box ------------------------------------------------
578 
579  xmin = minval(pgrid%xyz(xcoord,1:pgrid%nVertTot))
580  xmax = maxval(pgrid%xyz(xcoord,1:pgrid%nVertTot))
581  ymin = minval(pgrid%xyz(ycoord,1:pgrid%nVertTot))
582  ymax = maxval(pgrid%xyz(ycoord,1:pgrid%nVertTot))
583  zmin = minval(pgrid%xyz(zcoord,1:pgrid%nVertTot))
584  zmax = maxval(pgrid%xyz(zcoord,1:pgrid%nVertTot))
585 
586  xdel = xmax - xmin
587  ydel = ymax - ymin
588  zdel = zmax - zmin
589 
590  xmin = xmin - delfrac*xdel
591  xmax = xmax + delfrac*xdel
592  ymin = ymin - delfrac*ydel
593  ymax = ymax + delfrac*ydel
594  zmin = zmin - delfrac*zdel
595  zmax = zmax + delfrac*zdel
596 
597 ! ----- Create and build Octree -----------------------------------------------
598 
599  CALL rflu_createoctree(global,pgrid%nCellsTot)
600  CALL rflu_buildoctree(pgrid%cofg(xcoord,1:pgrid%nCellsTot), &
601  pgrid%cofg(ycoord,1:pgrid%nCellsTot), &
602  pgrid%cofg(zcoord,1:pgrid%nCellsTot), &
603  xmin,xmax,ymin,ymax,zmin,zmax)
604 
605 ! ----- Loop over faces -------------------------------------------------------
606 
607  DO ifc = 1,pgrid%nFaces
608  CALL rflu_queryoctree(pgrid%fc(xcoord,ifc),pgrid%fc(ycoord,ifc), &
609  pgrid%fc(zcoord,ifc),stencilsizetemp, &
610  stenciltemp)
611 
612 ! ------- Compute distance (signed) normal to face. At present, this only
613 ! accepts cells if they are located in the normal direction of the face,
614 ! so you can only generate 1x1x<nCells> stencils with this method. Not
615 ! hard to get other stencils, but hard to control them on non-uniform
616 ! grids. Rejected faces get a distance value of HUGE(1.0_RFREAL), which
617 ! is used further below in checking that enough cells were found.
618 
619  dummy = maxloc(abs(pgrid%fn(1:3,ifc)))
620  loc = dummy(1)
621 
622  DO is = 1,stencilsizetemp
623  ic = stenciltemp(is)
624 
625  vec(:) = pgrid%cofg(:,ic) - pgrid%fc(:,ifc)
626  vecm = sqrt(vec(1)*vec(1) + vec(2)*vec(2) + vec(3)*vec(3))
627  vec(:) = vec(:)/vecm
628 
629  dotp = vec(1)*pgrid%fn(1,ifc) &
630  + vec(2)*pgrid%fn(2,ifc) &
631  + vec(3)*pgrid%fn(3,ifc)
632 
633  IF ( floatequal(abs(dotp),1.0_rfreal) .EQV. .true. ) THEN
634  dist(is) = pgrid%cofg(loc,ic) - pgrid%fc(loc,ifc)
635  ELSE
636  dist(is) = huge(1.0_rfreal)
637  END IF ! dist
638  END DO ! is
639 
640 ! ------- Sort cells according to increasing distance
641 
642  CALL quicksortrfrealinteger(dist(1:stencilsizetemp), &
643  stenciltemp(1:stencilsizetemp), &
644  stencilsizetemp)
645 
646 ! ------- Make sure that stencil contains cells which straddle face, and that
647 ! the two cells are in the correct location (very strict check).
648 
649  IF ( pgrid%fn(loc,ifc) > 0.0_rfreal ) THEN
650  locoffset1 = 0
651  locoffset2 = 1
652  ELSE
653  locoffset1 = 1
654  locoffset2 = 0
655  END IF ! pGrid%fn
656 
657  DO ic = 1,2
658  IF ( ic == 1 ) THEN ! NOTE integer division
659  loc = stencilsize/2 + locoffset1
660  ELSE
661  loc = stencilsize/2 + locoffset2
662  END IF ! ic
663 
664  IF ( stenciltemp(loc) /= pgrid%f2c(ic,ifc) ) THEN
665  CALL errorstop(global,err_oles_stencil,__line__)
666  END IF ! stencilTemp
667  END DO ! ic
668 
669 ! ------- Make sure that the first <stencilSize> cells are in line with the face,
670 ! in other words, do not have dist(is) = HUGE(1.0_RFREAL)
671 
672  DO is = 1,stencilsize
673  IF ( floatequal(dist(is),huge(1.0_rfreal)) .EQV. .true. ) THEN
674  CALL errorstop(global,err_oles_stencil,__line__)
675  END IF ! dist
676  END DO ! is
677 
678 ! ------- Copy first <stencilSize> cells into face stencil array
679 
680  DO is = 1,stencilsize
681  pgrid%fsOLES(is,ifc) = stenciltemp(is)
682  END DO ! is
683  END DO ! ifc
684 
685 ! ----- Destroy Octree and deallocate memory ----------------------------------------
686 
687  CALL rflu_destroyoctree(global)
688 
689  DEALLOCATE(stenciltemp,stat=errorflag)
690  global%error = errorflag
691  IF ( global%error /= err_none ) THEN
692  CALL errorstop(global,err_deallocate,__line__,'stencilTemp')
693  END IF ! global%error
694 
695  DEALLOCATE(dist,stat=errorflag)
696  global%error = errorflag
697  IF ( global%error /= err_none ) THEN
698  CALL errorstop(global,err_deallocate,__line__,'dist')
699  END IF ! global%error
700  END IF ! stencilSize
701 
702 
703 #ifdef CHECK_DATASTRUCT
704 ! --- Data structure output for checking
705  WRITE(stdout,'(A)') solver_name
706  WRITE(stdout,'(A,1X,A)') solver_name,'### START CHECK OUTPUT ###'
707  WRITE(stdout,'(A,1X,A)') solver_name,'Optimal LES face stencils'
708  IF ( stencilsize == 2 ) THEN
709  DO ifc = 1,pgrid%nFaces
710  WRITE(stdout,'(A,3(1X,I6))') solver_name,ifc,pgrid%fsOLES(1:2,ifc)
711  END DO ! ifc
712  ELSE IF ( stencilsize == 4 ) THEN
713  DO ifc = 1,pgrid%nFaces
714  WRITE(stdout,'(A,5(1X,I6))') solver_name,ifc,pgrid%fsOLES(1:4,ifc)
715  END DO ! ifc
716  END IF ! stencilSize
717  WRITE(stdout,'(A,1X,A)') solver_name,'### END CHECK OUTPUT ###'
718  WRITE(stdout,'(A)') solver_name
719 #endif
720 
721 ! ==============================================================================
722 ! End
723 ! ==============================================================================
724 
725  IF ( global%myProcid == masterproc .AND. &
726  global%verbLevel > verbose_none ) THEN
727  WRITE(stdout,'(A,1X,A)') solver_name,'Building optimal LES '// &
728  'stencils done.'
729  END IF ! global%verbLevel
730 
731  CALL deregisterfunction(global)
732 
733  END SUBROUTINE rflu_buildstencilsoles
734 
735 
736 
737 
738 ! ******************************************************************************
739 ! Find prototypical faces: to reduce computational cost
740 ! ******************************************************************************
741 
742  SUBROUTINE rflu_findprototypefacesoles(pRegion)
743 
744  IMPLICIT NONE
745 
746 ! --- parameters
747 
748  TYPE(t_region), POINTER :: pregion
749 
750 ! --- local variables
751 
752  INTEGER, PARAMETER :: locsave_init = -1 ! Must be zero or negative
753  INTEGER :: ifc,ifcp,loc,locsavecntr
754  INTEGER :: dummy(1),locsave(3)
755  TYPE(t_grid), POINTER :: pgrid
756  TYPE(t_global), POINTER :: global
757 
758 ! ==============================================================================
759 ! Start
760 ! ==============================================================================
761 
762  global => pregion%global
763 
764  CALL registerfunction(global,'RFLU_FindPrototypeFacesOLES',&
765  'RFLU_ModOLES.F90')
766 
767  IF ( global%myProcid == masterproc .AND. &
768  global%verbLevel > verbose_none ) THEN
769  WRITE(stdout,'(A)') solver_name
770  WRITE(stdout,'(A,1X,A)') solver_name,'Finding optimal LES '// &
771  'prototype faces...'
772  END IF ! global%verbLevel
773 
774 ! ==============================================================================
775 ! Set grid pointer
776 ! ==============================================================================
777 
778  pgrid => pregion%grid
779 
780 ! ==============================================================================
781 ! Initialize
782 ! ==============================================================================
783 
784  locsavecntr = 0
785  locsave(1:3) = locsave_init
786 
787 ! ==============================================================================
788 ! Find prototype faces
789 ! ==============================================================================
790 
791  IF ( global%myProcid == masterproc .AND. &
792  global%verbLevel > verbose_none ) THEN
793  WRITE(stdout,'(A,3X,A)') solver_name,'Finding prototype faces...'
794  END IF ! global%verbLevel
795 
796  DO ifc = 1,pgrid%nFaces
797  dummy = maxloc(abs(pgrid%fn(1:3,ifc)))
798  loc = dummy(1)
799 
800  IF ( locsave(loc) == locsave_init ) THEN
801  locsave(loc) = 1
802  locsavecntr = locsavecntr + 1
803  pgrid%fp2fOLES(loc) = ifc
804  END IF ! locSave
805 
806  IF ( locsavecntr == 3 ) THEN
807  EXIT
808  END IF ! locSaveCntr
809  END DO ! ifc
810 
811 #ifdef CHECK_DATASTRUCT
812 ! --- Data structure output for checking
813  WRITE(stdout,'(A)') solver_name
814  WRITE(stdout,'(A,1X,A)') solver_name,'### START CHECK OUTPUT ###'
815  WRITE(stdout,'(A,1X,A)') solver_name,'Prototype faces:'
816  DO ifcp = 1,3
817  ifc = pgrid%fp2fOLES(ifcp)
818  WRITE(stdout,'(A,2(1X,I6),3(1X,E18.9))') solver_name,ifcp,ifc, &
819  pgrid%fn(1:3,ifc)
820  END DO ! ifcp
821  WRITE(stdout,'(A,1X,A)') solver_name,'### END CHECK OUTPUT ###'
822  WRITE(stdout,'(A)') solver_name
823 #endif
824 
825 ! ==============================================================================
826 ! Map other faces onto prototype faces
827 ! ==============================================================================
828 
829  IF ( global%myProcid == masterproc .AND. &
830  global%verbLevel > verbose_none ) THEN
831  WRITE(stdout,'(A,3X,A)') solver_name,'Mapping faces to '// &
832  'prototype faces...'
833  END IF ! global%verbLevel
834 
835  DO ifc = 1,pgrid%nFaces
836  dummy = maxloc(abs(pgrid%fn(1:3,ifc)))
837  loc = dummy(1)
838 
839  pgrid%f2fpOLES(ifc) = loc
840  END DO ! ifc
841 
842 #ifdef CHECK_DATASTRUCT
843 ! --- Data structure output for checking
844  WRITE(stdout,'(A)') solver_name
845  WRITE(stdout,'(A,1X,A)') solver_name,'### START CHECK OUTPUT ###'
846  WRITE(stdout,'(A,1X,A)') solver_name,'Prototype faces:'
847  DO ifc = 1,pgrid%nFaces
848  ifcp = pgrid%f2fpOLES(ifc)
849  WRITE(stdout,'(A,2(1X,I6),3(1X,E18.9))') solver_name,ifc,ifcp, &
850  pgrid%fn(1:3,ifc)
851  END DO ! ifc
852  WRITE(stdout,'(A,1X,A)') solver_name,'### END CHECK OUTPUT ###'
853  WRITE(stdout,'(A)') solver_name
854 #endif
855 
856 ! ==============================================================================
857 ! End
858 ! ==============================================================================
859 
860  IF ( global%myProcid == masterproc .AND. &
861  global%verbLevel > verbose_none ) THEN
862  WRITE(stdout,'(A,1X,A)') solver_name,'Finding optimal LES '// &
863  'prototype faces done.'
864  END IF ! global%verbLevel
865 
866  CALL deregisterfunction(global)
867 
868  END SUBROUTINE rflu_findprototypefacesoles
869 
870 
871 
872 
873 
874 
875 
876 
877 ! ******************************************************************************
878 ! Build stencil
879 ! ******************************************************************************
880 
881  SUBROUTINE rflu_computegeometrictermsoles(pRegion)
882 
883  IMPLICIT NONE
884 
885 ! --- parameters
886 
887  TYPE(t_region), POINTER :: pregion
888 
889 ! --- local variables
890 
891  INTEGER :: ic,icl,icg,ifc,ip,c1,c2,loc,ncells
892  INTEGER :: dummy(1)
893  REAL(RFREAL) :: term
894  REAL(RFREAL) :: fc(xcoord:zcoord)
895  REAL(RFREAL) :: xyz(int_lim_low:int_lim_upp,xcoord:zcoord)
896  TYPE(t_grid), POINTER :: pgrid
897  TYPE(t_patch), POINTER :: ppatch
898  TYPE(t_global), POINTER :: global
899 
900 #ifdef CHECK_DATASTRUCT
901  INTEGER :: i,j
902 #endif
903 
904 ! ==============================================================================
905 ! Start
906 ! ==============================================================================
907 
908  global => pregion%global
909 
910  CALL registerfunction(global,'RFLU_ComputeGeometricTermsOLES',&
911  'RFLU_ModOLES.F90')
912 
913  IF ( global%myProcid == masterproc .AND. &
914  global%verbLevel > verbose_none ) THEN
915  WRITE(stdout,'(A)') solver_name
916  WRITE(stdout,'(A,1X,A)') solver_name,'Computing optimal LES '// &
917  'geometric terms...'
918  END IF ! global%verbLevel
919 
920 ! ==============================================================================
921 ! Set grid pointer
922 ! ==============================================================================
923 
924  pgrid => pregion%grid
925 
926 ! ==============================================================================
927 ! Compute integration limits for each cell
928 ! ==============================================================================
929 
930  IF ( global%myProcid == masterproc .AND. &
931  global%verbLevel > verbose_none ) THEN
932  WRITE(stdout,'(A,3X,A)') solver_name,'Computing cell integration '// &
933  'limits...'
934  IF ( global%verbLevel > verbose_low ) THEN
935  WRITE(stdout,'(A,5X,A)') solver_name,'Interior faces...'
936  END IF ! global%verbLevel
937  END IF ! global%verbLevel
938 
939  DO ifc = 1,pgrid%nFacesTot ! NOTE loop over ALL faces
940  c1 = pgrid%f2c(1,ifc)
941  c2 = pgrid%f2c(2,ifc)
942 
943  dummy = maxloc(abs(pgrid%fn(1:3,ifc)))
944  loc = dummy(1)
945 
946  fc(:) = abs(pgrid%fn(loc,ifc))*pgrid%fc(:,ifc)
947 
948  IF ( c2 /= f2c_init ) THEN
949  IF ( pgrid%fn(loc,ifc) > 0.0_rfreal ) THEN
950  pgrid%intLimOLES(int_lim_upp,loc,c1) = fc(loc)
951  pgrid%intLimOLES(int_lim_low,loc,c2) = fc(loc)
952  ELSE
953  pgrid%intLimOLES(int_lim_low,loc,c1) = fc(loc)
954  pgrid%intLimOLES(int_lim_upp,loc,c2) = fc(loc)
955  END IF ! pGrid%intLimOLES
956  ELSE
957  IF ( pgrid%fn(loc,ifc) > 0.0_rfreal ) THEN
958  pgrid%intLimOLES(int_lim_upp,loc,c1) = fc(loc)
959  ELSE
960  pgrid%intLimOLES(int_lim_low,loc,c1) = fc(loc)
961  END IF ! pGrid%intLimOLES
962  END IF ! c2
963  END DO ! ifc
964 
965  IF ( global%myProcid == masterproc .AND. &
966  global%verbLevel > verbose_low ) THEN
967  WRITE(stdout,'(A,5X,A)') solver_name,'Patches...'
968  END IF ! global%verbLevel
969 
970  DO ip = 1,pgrid%nPatches
971  ppatch => pregion%patches(ip)
972 
973  IF ( global%myProcid == masterproc .AND. &
974  global%verbLevel > verbose_low ) THEN
975  WRITE(stdout,'(A,7X,A,I4)') solver_name,'Patch: ',ip
976  END IF ! global%verbLevel
977 
978  DO ifc = 1,ppatch%nBFaces
979  c1 = ppatch%bf2c(ifc)
980 
981  dummy = maxloc(abs(ppatch%fn(1:3,ifc)))
982  loc = dummy(1)
983 
984  fc(:) = abs(ppatch%fn(loc,ifc))*ppatch%fc(:,ifc)
985 
986  IF ( ppatch%fn(loc,ifc) > 0.0_rfreal ) THEN
987  pgrid%intLimOLES(int_lim_upp,loc,c1) = fc(loc)
988  ELSE
989  pgrid%intLimOLES(int_lim_low,loc,c1) = fc(loc)
990  END IF ! pGrid%intLimOLES
991  END DO ! ifc
992  END DO ! ip
993 
994 
995 #ifdef CHECK_DATASTRUCT
996 ! --- Data structure output for checking
997  WRITE(stdout,'(A)') solver_name
998  WRITE(stdout,'(A,1X,A)') solver_name,'### START CHECK OUTPUT ###'
999  WRITE(stdout,'(A,1X,A)') solver_name,'Optimal LES cell integration limits'
1000  DO ic = 1,pgrid%nCellsTot
1001  WRITE(stdout,'(A,1X,I6,6(1X,E18.9))') solver_name,ic, &
1002  ((pgrid%intLimOLES(i,j,ic),i=int_lim_low,int_lim_upp),j=xcoord,zcoord)
1003  END DO ! ic
1004  WRITE(stdout,'(A,1X,A)') solver_name,'### END CHECK OUTPUT ###'
1005  WRITE(stdout,'(A)') solver_name
1006 #endif
1007 
1008 ! ==============================================================================
1009 ! Compute radius of smallest sphere contained in cells of stencil
1010 ! ==============================================================================
1011 
1012  IF ( global%myProcid == masterproc .AND. &
1013  global%verbLevel > verbose_none ) THEN
1014  WRITE(stdout,'(A,3X,A)') solver_name,'Computing smallest sphere '// &
1015  'for stencils...'
1016  END IF ! global%verbLevel
1017 
1018  DO ifc = 1,pgrid%nFaces
1019  ncells = SIZE(pgrid%fsOLES,1)
1020 
1021  xyz(int_lim_low,xcoord:zcoord) = huge(1.0_rfreal)
1022  xyz(int_lim_upp,xcoord:zcoord) = -huge(1.0_rfreal)
1023 
1024  DO icl = 1,ncells
1025  icg = pgrid%fsOLES(icl,ifc)
1026 
1027  xyz(int_lim_low,1:3) = min(pgrid%intLimOLES(int_lim_low,1:3,icg), &
1028  xyz(int_lim_low,1:3))
1029  xyz(int_lim_upp,1:3) = max(pgrid%intLimOLES(int_lim_upp,1:3,icg), &
1030  xyz(int_lim_upp,1:3))
1031  END DO ! ic
1032 
1033 ! BEGIN ALTERNATIVE 1: Wrong, represents sphere containED in cells
1034  DO ic = 1,3
1035  pgrid%rhoOLES(ifc) = min(pgrid%rhoOLES(ifc), &
1036  xyz(int_lim_upp,ic)-xyz(int_lim_low,ic))
1037  END DO ! ic
1038 ! END ALTERNATIVE 1
1039 
1040 ! BEGIN ALERNATIVE 2: Correct, represents sphere containing cells
1041 ! Changes integrals, but not weights
1042 ! term = -HUGE(1.0_RFREAL)
1043 !
1044 ! DO ic = 1,3
1045 ! term = MAX(term,xyz(INT_LIM_UPP,ic)-xyz(INT_LIM_LOW,ic))
1046 ! END DO ! ic
1047 !
1048 ! pGrid%rhoOLES(ifc) = MIN(pGrid%rhoOLES(ifc),term)
1049 ! END ALTERNATIVE 2
1050  END DO ! ifc
1051 
1052 #ifdef CHECK_DATASTRUCT
1053 ! --- Data structure output for checking
1054  WRITE(stdout,'(A)') solver_name
1055  WRITE(stdout,'(A,1X,A)') solver_name,'### START CHECK OUTPUT ###'
1056  WRITE(stdout,'(A,1X,A)') solver_name,'Optimal LES smallest stencil sphere'
1057  DO ifc = 1,pgrid%nFaces
1058  WRITE(stdout,'(A,1X,I6,1X,E18.9)') solver_name,ifc,pgrid%rhoOLES(ifc)
1059  END DO ! ifc
1060  WRITE(stdout,'(A,1X,A)') solver_name,'### END CHECK OUTPUT ###'
1061  WRITE(stdout,'(A)') solver_name
1062 #endif
1063 
1064 ! ==============================================================================
1065 ! Compute filter width - assume constant grid spacing for now
1066 ! ==============================================================================
1067 
1068  IF ( global%myProcid == masterproc .AND. &
1069  global%verbLevel > verbose_none ) THEN
1070  WRITE(stdout,'(A,3X,A)') solver_name,'Computing filter width...'
1071  END IF ! global%verbLevel
1072 
1073  pgrid%deltaOLES = pgrid%vol(1)**(1.0_rfreal/3.0_rfreal)
1074 
1075 #ifdef CHECK_DATASTRUCT
1076 ! --- Data structure output for checking
1077  WRITE(stdout,'(A)') solver_name
1078  WRITE(stdout,'(A,1X,A)') solver_name,'### START CHECK OUTPUT ###'
1079  WRITE(stdout,'(A,1X,A)') solver_name,'Optimal LES filter width'
1080  WRITE(stdout,'(A,1X,E18.9)') solver_name,pgrid%deltaOLES
1081  WRITE(stdout,'(A,1X,A)') solver_name,'### END CHECK OUTPUT ###'
1082  WRITE(stdout,'(A)') solver_name
1083 #endif
1084 
1085 ! ==============================================================================
1086 ! End
1087 ! ==============================================================================
1088 
1089  IF ( global%myProcid == masterproc .AND. &
1090  global%verbLevel > verbose_none ) THEN
1091  WRITE(stdout,'(A,1X,A)') solver_name,'Computing optimal LES '// &
1092  'geometric terms done.'
1093  END IF ! global%verbLevel
1094 
1095  CALL deregisterfunction(global)
1096 
1097  END SUBROUTINE rflu_computegeometrictermsoles
1098 
1099 
1100 
1101 
1102 ! ******************************************************************************
1103 ! Build symmetry maps for integrals
1104 ! ******************************************************************************
1105 
1106  SUBROUTINE rflu_buildsymmetrymapsoles(pRegion)
1107 
1108  IMPLICIT NONE
1109 
1110 ! --- parameters
1111 
1112  TYPE(t_region), POINTER :: pregion
1113 
1114 ! --- local variables
1115 
1116  INTEGER :: b,g,j,k,ncells,pos,possym
1117  TYPE(t_global), POINTER :: global
1118 
1119 ! ==============================================================================
1120 ! Start
1121 ! ==============================================================================
1122 
1123  global => pregion%global
1124 
1125  CALL registerfunction(global,'RFLU_BuildSymmetryMapsOLES',&
1126  'RFLU_ModOLES.F90')
1127 
1128  IF ( global%myProcid == masterproc .AND. &
1129  global%verbLevel > verbose_none ) THEN
1130  WRITE(stdout,'(A)') solver_name
1131  WRITE(stdout,'(A,1X,A)') solver_name,'Building optimal LES '// &
1132  'symmetry maps...'
1133  END IF ! global%verbLevel
1134 
1135 ! ==============================================================================
1136 ! Build symmetry maps
1137 ! ==============================================================================
1138 
1139  ncells = SIZE(pregion%grid%fsOLES,1)
1140 
1141  DO b = 1,ncells
1142  DO g = 1,ncells
1143  DO j = 1,3
1144  DO k = 1,3
1145  pos = rflu_getqposoles(j,k,b,g,ncells)
1146  possym = rflu_getqposoles(k,j,g,b,ncells)
1147 
1148  IF ( possym < pos ) THEN
1149  symmapi45oles(pos) = possym
1150  ELSE
1151  symmapi45oles(pos) = 0
1152  END IF ! pos2
1153  END DO ! k
1154  END DO ! j
1155  END DO ! g
1156  END DO ! b
1157 
1158 ! ==============================================================================
1159 ! End
1160 ! ==============================================================================
1161 
1162  IF ( global%myProcid == masterproc .AND. &
1163  global%verbLevel > verbose_none ) THEN
1164  WRITE(stdout,'(A,1X,A)') solver_name,'Building optimal LES '// &
1165  'symmetry maps done.'
1166  END IF ! global%verbLevel
1167 
1168  CALL deregisterfunction(global)
1169 
1170  END SUBROUTINE rflu_buildsymmetrymapsoles
1171 
1172 
1173 
1174 
1175 ! ******************************************************************************
1176 ! Enforce symmetry: NOTE will not be completely symmetric!
1177 ! ******************************************************************************
1178 
1179  SUBROUTINE rflu_enforcesymmetryoles(pRegion)
1180 
1181  IMPLICIT NONE
1182 
1183 ! --- parameters
1184 
1185  TYPE(t_region), POINTER :: pregion
1186 
1187 ! --- local variables
1188 
1189  INTEGER :: dloc,ifcp,hloc,ncells,vloc
1190  TYPE(t_global), POINTER :: global
1191 
1192 ! ==============================================================================
1193 ! Start
1194 ! ==============================================================================
1195 
1196  global => pregion%global
1197 
1198  CALL registerfunction(global,'RFLU_EnforceSymmetryOLES',&
1199  'RFLU_ModOLES.F90')
1200 
1201  IF ( global%myProcid == masterproc .AND. &
1202  global%verbLevel > verbose_none ) THEN
1203  WRITE(stdout,'(A)') solver_name
1204  WRITE(stdout,'(A,1X,A)') solver_name,'Enforcing optimal LES '// &
1205  'integral symmetry...'
1206  END IF ! global%verbLevel
1207 
1208 ! ==============================================================================
1209 ! Set grid pointer and variables
1210 ! ==============================================================================
1211 
1212  pgrid => pregion%grid
1213 
1214  ncells = SIZE(pgrid%fsOLES,1)
1215 
1216 ! ==============================================================================
1217 ! Enforce symmetry
1218 ! ==============================================================================
1219 
1220 ! ------------------------------------------------------------------------------
1221 ! Integral I4
1222 ! ------------------------------------------------------------------------------
1223 
1224  IF ( global%myProcid == masterproc .AND. &
1225  global%verbLevel > verbose_none ) THEN
1226  WRITE(stdout,'(A,3X,A)') solver_name,'Integral 4...'
1227  END IF ! global%verbLevel
1228 
1229  DO ifcp = 1,3
1230  DO vloc = 1,9*ncells*ncells
1231  IF ( symmapi45oles(vloc) /= 0 ) THEN
1232  pgrid%int40OLES(:,ifcp,vloc) = &
1233  pgrid%int40OLES(:,ifcp,symmapi45oles(vloc))
1234  pgrid%int41OLES(:,ifcp,vloc) = &
1235  pgrid%int41OLES(:,ifcp,symmapi45oles(vloc))
1236  pgrid%int42OLES(:,ifcp,vloc) = &
1237  pgrid%int42OLES(:,ifcp,symmapi45oles(vloc))
1238  END IF ! symMapI45OLES
1239  END DO ! vLoc
1240  END DO ! ifcp
1241 
1242 ! ------------------------------------------------------------------------------
1243 ! Integral I5
1244 ! ------------------------------------------------------------------------------
1245 
1246  IF ( global%myProcid == masterproc .AND. &
1247  global%verbLevel > verbose_none ) THEN
1248  WRITE(stdout,'(A,3X,A)') solver_name,'Integral 5...'
1249  END IF ! global%verbLevel
1250 
1251  DO ifcp = 1,3 ! equality of rows
1252  DO vloc = 1,9*ncells*ncells
1253  IF ( symmapi45oles(vloc) /= 0 ) THEN
1254  pgrid%int50OLES(ifcp,vloc,:) = &
1255  pgrid%int50OLES(ifcp,symmapi45oles(vloc),:)
1256  pgrid%int51OLES(ifcp,vloc,:) = &
1257  pgrid%int51OLES(ifcp,symmapi45oles(vloc),:)
1258  pgrid%int52OLES(ifcp,vloc,:) = &
1259  pgrid%int52OLES(ifcp,symmapi45oles(vloc),:)
1260  END IF ! symMapI45OLES
1261  END DO ! vLoc
1262  END DO ! ifcp
1263 
1264 ! This is screwing up the matrix - no longer singular (yes, should be because of
1265 ! symmetry of weights). This is strange, but have not yet found out why destroying
1266 ! singularity. Enforcing the symmetries should make it as singular as possible...
1267 ! DO ifcp = 1,3 ! symmetry of entire matrix
1268 ! DO dLoc = 1,9*nCells*nCells
1269 ! DO vLoc = dLoc+1,9*nCells*nCells
1270 ! hLoc = vLoc
1271 !
1272 ! pGrid%int50OLES(ifcp,dLoc,vLoc) = pGrid%int50OLES(ifcp,vLoc,dLoc)
1273 ! pGrid%int51OLES(ifcp,dLoc,vLoc) = pGrid%int51OLES(ifcp,vLoc,dLoc)
1274 ! pGrid%int52OLES(ifcp,dLoc,vLoc) = pGrid%int52OLES(ifcp,vLoc,dLoc)
1275 ! END DO ! hLoc
1276 ! END DO ! dLoc
1277 ! END DO ! ifcp
1278 
1279 ! ==============================================================================
1280 ! End
1281 ! ==============================================================================
1282 
1283  IF ( global%myProcid == masterproc .AND. &
1284  global%verbLevel > verbose_none ) THEN
1285  WRITE(stdout,'(A,1X,A)') solver_name,'Enforcing optimal LES '// &
1286  'integral symmetry done.'
1287  END IF ! global%verbLevel
1288 
1289  CALL deregisterfunction(global)
1290 
1291  END SUBROUTINE rflu_enforcesymmetryoles
1292 
1293 
1294 
1295 
1296 ! ******************************************************************************
1297 ! Define the correlation functions which are to be integrated over the cells
1298 ! ******************************************************************************
1299 
1300 ! ==============================================================================
1301 ! Second-order two-point correlation (for I^2), part 0
1302 ! ==============================================================================
1303 
1304  SUBROUTINE rflu_definecorrelation220(nDim,z,nFunNZ,f)
1305 
1306 ! --- Arguments
1307 
1308  INTEGER, INTENT(IN) :: ndim,nfunnz
1309  DOUBLE PRECISION, INTENT(IN) :: z(ndim)
1310  DOUBLE PRECISION, INTENT(INOUT) :: f(nfunnz)
1311 
1312 ! --- Locals
1313 
1314  INTEGER :: ifun,ifunnz,j,l
1315 
1316 ! --- Start
1317 
1318  DO ifunnz = 1,nfunnz
1319  ifun = mapfunnz2funcorr22(ifunnz)
1320  CALL rflu_mapk2ij(ifun,l,j)
1321 
1322  f(ifunnz) = kd(l,j)
1323  END DO ! iFunNZ
1324 
1325 ! --- End
1326 
1327  END SUBROUTINE rflu_definecorrelation220
1328 
1329 ! ==============================================================================
1330 ! Second-order two-point correlation (for I^2), part 1
1331 ! ==============================================================================
1332 
1333  SUBROUTINE rflu_definecorrelation221(nDim,z,nFunNZ,f)
1334 
1335 ! --- Arguments
1336 
1337  INTEGER, INTENT(IN) :: ndim,nfunnz
1338  DOUBLE PRECISION, INTENT(IN) :: z(ndim)
1339  DOUBLE PRECISION, INTENT(INOUT) :: f(nfunnz)
1340 
1341 ! --- Locals
1342 
1343  INTEGER :: ifun,ifunnz,j,l
1344  DOUBLE PRECISION :: fact,rm2,rm23,rm2i
1345 
1346 ! --- Start
1347 
1348  rm2 = (z(4)-z(1))**2 + (z(5)-z(2))**2 + (z(6)-z(3))**2
1349  rm23 = rm2**(1.0d0/3.0d0)
1350  rm2i = 1.0_rfreal/makenonzero(rm2)
1351 
1352  DO ifunnz = 1,nfunnz
1353  ifun = mapfunnz2funcorr22(ifunnz)
1354  CALL rflu_mapk2ij(ifun,l,j)
1355 
1356  f(ifunnz) = rm23*(rm2i*(z(l+3)-z(l))*(z(j+3)-z(j)) - 4.0d0*kd(l,j))
1357  END DO ! iFunNZ
1358 
1359 ! --- End
1360 
1361  END SUBROUTINE rflu_definecorrelation221
1362 
1363 
1364 
1365 ! ==============================================================================
1366 ! Third-order two-point correlation (for I^1)
1367 ! ==============================================================================
1368 
1369  SUBROUTINE rflu_definecorrelation32(nDim,z,nFunNZ,f)
1370 
1371 ! --- Arguments
1372 
1373  INTEGER, INTENT(IN) :: ndim,nfunnz
1374  DOUBLE PRECISION, INTENT(IN) :: z(ndim)
1375  DOUBLE PRECISION, INTENT(INOUT) :: f(nfunnz)
1376 
1377 ! --- Locals
1378 
1379  INTEGER :: i,ifun,ifunnz,l,n
1380  DOUBLE PRECISION :: term1,term2,term3
1381  DOUBLE PRECISION :: zd(ndim+1)
1382 
1383 ! --- Start
1384 
1385  n = nzloc
1386 
1387 ! --- Copy integration points and add constant face coordinate
1388 
1389  zd(1:ndim) = z(1:ndim)
1390  zd(ndim+1) = nzval
1391 
1392 ! --- Evaluate integral
1393 
1394  DO ifunnz = 1,nfunnz
1395  ifun = mapfunnz2funcorr32(ifunnz)
1396  CALL rflu_mapk2ij(ifun,l,i)
1397 
1398  term1 = kd(i,n)*(zd(l)-zd(mapsurf2vol2(n,l)))
1399  term2 = kd(i,l)*(zd(n)-zd(mapsurf2vol2(n,n)))
1400  term3 = kd(n,l)*(zd(i)-zd(mapsurf2vol2(n,i)))
1401 
1402  f(ifunnz) = nzsgn*(term1 - 1.5d0*(term2 + term3))
1403  END DO ! iFunNZ
1404 
1405 ! --- End
1406 
1407  END SUBROUTINE rflu_definecorrelation32
1408 
1409 
1410 
1411 ! ==============================================================================
1412 ! Fourth-order three-point correlation (for I^4), part 0
1413 ! ==============================================================================
1414 
1415  SUBROUTINE rflu_definecorrelation430(nDim,z,nFunNZ,f)
1416 
1417 ! --- Arguments
1418 
1419  INTEGER, INTENT(IN) :: ndim,nfunnz
1420  DOUBLE PRECISION, INTENT(IN) :: z(ndim)
1421  DOUBLE PRECISION, INTENT(INOUT) :: f(nfunnz)
1422 
1423 ! --- Locals
1424 
1425  INTEGER :: i,ifun,ifunnz,l,m,n
1426  DOUBLE PRECISION :: term1,term2,term3
1427 
1428 ! --- Start
1429 
1430  n = nzloc
1431 
1432 ! --- Evaluate integral
1433 
1434  DO ifunnz = 1,nfunnz
1435  ifun = mapfunnz2funcorr43(ifunnz)
1436  CALL rflu_mapl2ijk(ifun,l,m,i)
1437 
1438  f(ifunnz) = nzsgn*(kd(l,m)*kd(i,n) + kd(l,i)*kd(m,n) + kd(l,n)*kd(m,i))
1439  END DO ! iFunNZ
1440 
1441 ! --- End
1442 
1443  END SUBROUTINE rflu_definecorrelation430
1444 
1445 
1446 ! ==============================================================================
1447 ! Fourth-order three-point correlation (for I^4), part 1
1448 ! ==============================================================================
1449 
1450  SUBROUTINE rflu_definecorrelation431(nDim,z,nFunNZ,f)
1451 
1452 ! --- Arguments
1453 
1454  INTEGER, INTENT(IN) :: ndim,nfunnz
1455  DOUBLE PRECISION, INTENT(IN) :: z(ndim)
1456  DOUBLE PRECISION, INTENT(INOUT) :: f(nfunnz)
1457 
1458 ! --- Locals
1459 
1460  INTEGER :: i,ifun,ifunnz,j,k,l,m
1461  DOUBLE PRECISION :: r1m2,r1m23,r1m2i,r2m2,r2m23,r2m2i,r3m2,r3m23,r3m2i, &
1462  r4m2,r4m23,r4m2i,r5m2,r5m23,r5m2i,term1,term2,term3, &
1463  term4,term5
1464  DOUBLE PRECISION :: zd(ndim+1)
1465 
1466 ! --- Start
1467 
1468  j = nzloc
1469 
1470 ! --- Copy integration points and add constant face coordinate
1471 
1472  zd(1:ndim) = z(1:ndim)
1473  zd(ndim+1) = nzval
1474 
1475 ! --- Evaluate integral
1476 
1477  r1m2 = (zd(1) - zd(4))**2 &
1478  + (zd(2) - zd(5))**2 &
1479  + (zd(3) - zd(6))**2
1480  r2m2 = (zd(1) - zd(mapsurf2vol3(j,1)))**2 &
1481  + (zd(2) - zd(mapsurf2vol3(j,2)))**2 &
1482  + (zd(3) - zd(mapsurf2vol3(j,3)))**2
1483  r3m2 = r2m2
1484  r4m2 = (zd(4) - zd(mapsurf2vol3(j,4)))**2 &
1485  + (zd(5) - zd(mapsurf2vol3(j,5)))**2 &
1486  + (zd(6) - zd(mapsurf2vol3(j,6)))**2
1487  r5m2 = r4m2
1488 
1489  r1m23 = r1m2**(1.0d0/3.0d0)
1490  r2m23 = r2m2**(1.0d0/3.0d0)
1491  r3m23 = r2m23
1492  r4m23 = r4m2**(1.0d0/3.0d0)
1493  r5m23 = r4m23
1494 
1495  r1m2i = 1.0d0/makenonzero(r1m2)
1496  r2m2i = 1.0d0/makenonzero(r2m2)
1497  r3m2i = r2m2i
1498  r4m2i = 1.0d0/makenonzero(r4m2)
1499  r5m2i = r4m2i
1500 
1501  DO ifunnz = 1,nfunnz
1502  ifun = mapfunnz2funcorr43(ifunnz)
1503  CALL rflu_mapl2ijk(ifun,l,m,i)
1504 
1505  term1 = kd(i,j)*r1m23*( r1m2i*(zd(l )-zd(l+3)) &
1506  *(zd(m )-zd(m+3)) &
1507  - 4.0d0*kd(l,m))
1508 
1509  term2 = kd(l,i)*r5m23*( r5m2i*(zd(m+3)-zd(mapsurf2vol3(j,m))) &
1510  *(zd(j+3)-zd(mapsurf2vol3(j,j))) &
1511  - 4.0d0*kd(m,j))
1512 
1513  term3 = kd(m,j)*r2m23*( r2m2i*(zd(l )-zd(mapsurf2vol3(j,l))) &
1514  *(zd(i )-zd(mapsurf2vol3(j,i))) &
1515  - 4.0d0*kd(l,i))
1516 
1517  term4 = kd(l,j)*r4m23*( r4m2i*(zd(m+3)-zd(mapsurf2vol3(j,m))) &
1518  *(zd(i+3)-zd(mapsurf2vol3(j,i))) &
1519  - 4.0d0*kd(m,i))
1520 
1521  term5 = kd(m,i)*r3m23*( r3m2i*(zd(l )-zd(mapsurf2vol3(j,l))) &
1522  *(zd(j )-zd(mapsurf2vol3(j,j))) &
1523  - 4.0d0*kd(l,j))
1524  f(ifunnz) = nzsgn*(term1 + term2 + term3 + term4 + term5)
1525  END DO ! iFunNZ
1526 
1527 ! --- End
1528 
1529  END SUBROUTINE rflu_definecorrelation431
1530 
1531 
1532 ! ==============================================================================
1533 ! Fourth-order three-point correlation (for I^4), part 2
1534 ! ==============================================================================
1535 
1536  SUBROUTINE rflu_definecorrelation432(nDim,z,nFunNZ,f)
1537 
1538 ! --- Arguments
1539 
1540  INTEGER, INTENT(IN) :: ndim,nfunnz
1541  DOUBLE PRECISION, INTENT(IN) :: z(ndim)
1542  DOUBLE PRECISION, INTENT(INOUT) :: f(nfunnz)
1543 
1544 ! --- Locals
1545 
1546  INTEGER :: i,ifun,ifunnz,j,k,l,m
1547  DOUBLE PRECISION :: r2m2,r2m23,r2m2i,r3m2,r3m23,r3m2i,r4m2,r4m23,r4m2i, &
1548  r5m2,r5m23,r5m2i,term1,term2
1549  DOUBLE PRECISION :: zd(ndim+1)
1550 
1551 ! --- Start
1552 
1553  j = nzloc
1554 
1555 ! --- Copy integration points and add constant face coordinate
1556 
1557  zd(1:ndim) = z(1:ndim)
1558  zd(ndim+1) = nzval
1559 
1560 ! --- Evaluate integral
1561 
1562  r2m2 = (zd(1) - zd(mapsurf2vol3(j,1)))**2 &
1563  + (zd(2) - zd(mapsurf2vol3(j,2)))**2 &
1564  + (zd(3) - zd(mapsurf2vol3(j,3)))**2
1565  r3m2 = r2m2
1566  r4m2 = (zd(4) - zd(mapsurf2vol3(j,4)))**2 &
1567  + (zd(5) - zd(mapsurf2vol3(j,5)))**2 &
1568  + (zd(6) - zd(mapsurf2vol3(j,6)))**2
1569  r5m2 = r4m2
1570 
1571  r2m23 = r2m2**(1.0d0/3.0d0)
1572  r3m23 = r2m23
1573  r4m23 = r4m2**(1.0d0/3.0d0)
1574  r5m23 = r4m23
1575 
1576  r2m2i = 1.0d0/makenonzero(r2m2)
1577  r3m2i = r2m2i
1578  r4m2i = 1.0d0/makenonzero(r4m2)
1579  r5m2i = r4m2i
1580 
1581  DO ifunnz = 1,nfunnz
1582  ifun = mapfunnz2funcorr43(ifunnz)
1583  CALL rflu_mapl2ijk(ifun,l,m,i)
1584 
1585  term1 = &
1586  r2m23*r5m23*(r2m2i*(zd(l ) - zd(mapsurf2vol3(j,l))) &
1587  *(zd(i ) - zd(mapsurf2vol3(j,i))) - 4.0d0*kd(l,i)) &
1588  *(r5m2i*(zd(m+3) - zd(mapsurf2vol3(j,m))) &
1589  *(zd(j+3) - zd(mapsurf2vol3(j,j))) - 4.0d0*kd(m,j))
1590 
1591  term2 = &
1592  r3m23*r4m23*(r3m2i*(zd(l ) - zd(mapsurf2vol3(j,l))) &
1593  *(zd(j ) - zd(mapsurf2vol3(j,j))) - 4.0d0*kd(l,j)) &
1594  *(r4m2i*(zd(m+3) - zd(mapsurf2vol3(j,m))) &
1595  *(zd(i+3) - zd(mapsurf2vol3(j,i))) - 4.0d0*kd(m,i))
1596 
1597  f(ifunnz) = nzsgn*(term1 + term2)
1598  END DO ! iFunNZ
1599 
1600 ! --- End
1601 
1602  END SUBROUTINE rflu_definecorrelation432
1603 
1604 
1605 
1606 
1607 
1608 ! ==============================================================================
1609 ! Fourth-order four-point correlation (for I^5), part 0
1610 ! ==============================================================================
1611 
1612  SUBROUTINE rflu_definecorrelation540(nDim,z,nFunNZ,f)
1613 
1614 ! --- Arguments
1615 
1616  INTEGER, INTENT(IN) :: ndim,nfunnz
1617  DOUBLE PRECISION, INTENT(IN) :: z(ndim)
1618  DOUBLE PRECISION, INTENT(INOUT) :: f(nfunnz)
1619 
1620 ! --- Locals
1621 
1622  INTEGER :: ifun,ifunnz,j,l,k,m
1623  DOUBLE PRECISION :: term1,term2,term3
1624 
1625 ! --- Start, evaluate integral
1626 
1627  DO ifunnz = 1,nfunnz
1628  ifun = mapfunnz2funcorr44(ifunnz)
1629  CALL rflu_mapm2ijkl(ifun,l,m,j,k)
1630 
1631  f(ifunnz) = kd(l,m)*kd(j,k) + kd(l,j)*kd(m,k) + kd(l,k)*kd(m,j)
1632  END DO ! iFunNZ
1633 
1634 ! --- End
1635 
1636  END SUBROUTINE rflu_definecorrelation540
1637 
1638 
1639 
1640 ! ==============================================================================
1641 ! Fourth-order four-point correlation (for I^5), part 1
1642 ! ==============================================================================
1643 
1644  SUBROUTINE rflu_definecorrelation541(nDim,z,nFunNZ,f)
1645 
1646 ! --- Arguments
1647 
1648  INTEGER, INTENT(IN) :: ndim,nfunnz
1649  DOUBLE PRECISION, INTENT(IN) :: z(ndim)
1650  DOUBLE PRECISION, INTENT(INOUT) :: f(nfunnz)
1651 
1652 ! --- Locals
1653 
1654  INTEGER :: ifun,ifunnz,j,k,l,m
1655  DOUBLE PRECISION :: r1m2,r1m23,r1m2i,r2m2,r2m23,r2m2i,r3m2,r3m23,r3m2i, &
1656  r4m2,r4m23,r4m2i,r5m2,r5m23,r5m2i,r6m2,r6m23,r6m2i, &
1657  term1,term2,term3,term4,term5,term6
1658 
1659 ! --- Start, evaluate integral
1660 
1661  r1m2 = (z(1) - z( 4))**2 &
1662  + (z(2) - z( 5))**2 &
1663  + (z(3) - z( 6))**2
1664  r2m2 = (z(1) - z( 7))**2 &
1665  + (z(2) - z( 8))**2 &
1666  + (z(3) - z( 9))**2
1667  r3m2 = (z(1) - z(10))**2 &
1668  + (z(2) - z(11))**2 &
1669  + (z(3) - z(12))**2
1670  r4m2 = (z(4) - z( 7))**2 &
1671  + (z(5) - z( 8))**2 &
1672  + (z(6) - z( 9))**2
1673  r5m2 = (z(4) - z(10))**2 &
1674  + (z(5) - z(11))**2 &
1675  + (z(6) - z(12))**2
1676  r6m2 = (z(7) - z(10))**2 &
1677  + (z(8) - z(11))**2 &
1678  + (z(9) - z(12))**2
1679 
1680  r1m23 = r1m2**(1.0d0/3.0d0)
1681  r2m23 = r2m2**(1.0d0/3.0d0)
1682  r3m23 = r3m2**(1.0d0/3.0d0)
1683  r4m23 = r4m2**(1.0d0/3.0d0)
1684  r5m23 = r5m2**(1.0d0/3.0d0)
1685  r6m23 = r6m2**(1.0d0/3.0d0)
1686 
1687  r1m2i = 1.0d0/makenonzero(r1m2)
1688  r2m2i = 1.0d0/makenonzero(r2m2)
1689  r3m2i = 1.0d0/makenonzero(r3m2)
1690  r4m2i = 1.0d0/makenonzero(r4m2)
1691  r5m2i = 1.0d0/makenonzero(r5m2)
1692  r6m2i = 1.0d0/makenonzero(r6m2)
1693 
1694  DO ifunnz = 1,nfunnz
1695  ifun = mapfunnz2funcorr44(ifunnz)
1696  CALL rflu_mapm2ijkl(ifun,l,m,j,k)
1697 
1698  term1 = r6m23*kd(l,m)*(r6m2i*(z(j+6) - z(j+9)) &
1699  *(z(k+6) - z(k+9)) - 4.0d0*kd(j,k))
1700  term2 = r1m23*kd(j,k)*(r1m2i*(z(l ) - z(l+3)) &
1701  *(z(m ) - z(m+3)) - 4.0d0*kd(l,m))
1702  term3 = r5m23*kd(l,j)*(r5m2i*(z(m+3) - z(m+9)) &
1703  *(z(k+3) - z(k+9)) - 4.0d0*kd(m,k))
1704  term4 = r2m23*kd(m,k)*(r2m2i*(z(l ) - z(l+6)) &
1705  *(z(j ) - z(j+6)) - 4.0d0*kd(l,j))
1706  term5 = r4m23*kd(l,k)*(r4m2i*(z(m+3) - z(m+6)) &
1707  *(z(j+3) - z(j+6)) - 4.0d0*kd(m,j))
1708  term6 = r3m23*kd(m,j)*(r3m2i*(z(l ) - z(l+9)) &
1709  *(z(k ) - z(k+9)) - 4.0d0*kd(l,k))
1710 
1711  f(ifunnz) = term1 + term2 + term3 + term4 + term5 + term6
1712  END DO ! iFunNZ
1713 
1714 ! --- End
1715 
1716  END SUBROUTINE rflu_definecorrelation541
1717 
1718 
1719 
1720 ! ==============================================================================
1721 ! Fourth-order four-point correlation (for I^5), part 2
1722 ! ==============================================================================
1723 
1724  SUBROUTINE rflu_definecorrelation542(nDim,z,nFunNZ,f)
1725 
1726 ! --- Arguments
1727 
1728  INTEGER, INTENT(IN) :: ndim,nfunnz
1729  DOUBLE PRECISION, INTENT(IN) :: z(ndim)
1730  DOUBLE PRECISION, INTENT(INOUT) :: f(nfunnz)
1731 
1732 ! --- Locals
1733 
1734  INTEGER :: ifun,ifunnz,j,k,l,m
1735  DOUBLE PRECISION :: r1m2,r1m23,r1m2i,r2m2,r2m23,r2m2i,r3m2,r3m23,r3m2i, &
1736  r4m2,r4m23,r4m2i,r5m2,r5m23,r5m2i,r6m2,r6m23,r6m2i, &
1737  term1,term2,term3
1738 
1739 ! --- Start, evaluate integral
1740 
1741  r1m2 = (z(1) - z( 4))**2 &
1742  + (z(2) - z( 5))**2 &
1743  + (z(3) - z( 6))**2
1744  r2m2 = (z(1) - z( 7))**2 &
1745  + (z(2) - z( 8))**2 &
1746  + (z(3) - z( 9))**2
1747  r3m2 = (z(1) - z(10))**2 &
1748  + (z(2) - z(11))**2 &
1749  + (z(3) - z(12))**2
1750  r4m2 = (z(4) - z( 7))**2 &
1751  + (z(5) - z( 8))**2 &
1752  + (z(6) - z( 9))**2
1753  r5m2 = (z(4) - z(10))**2 &
1754  + (z(5) - z(11))**2 &
1755  + (z(6) - z(12))**2
1756  r6m2 = (z(7) - z(10))**2 &
1757  + (z(8) - z(11))**2 &
1758  + (z(9) - z(12))**2
1759 
1760  r1m23 = r1m2**(1.0d0/3.0d0)
1761  r2m23 = r2m2**(1.0d0/3.0d0)
1762  r3m23 = r3m2**(1.0d0/3.0d0)
1763  r4m23 = r4m2**(1.0d0/3.0d0)
1764  r5m23 = r5m2**(1.0d0/3.0d0)
1765  r6m23 = r6m2**(1.0d0/3.0d0)
1766 
1767  r1m2i = 1.0d0/makenonzero(r1m2)
1768  r2m2i = 1.0d0/makenonzero(r2m2)
1769  r3m2i = 1.0d0/makenonzero(r3m2)
1770  r4m2i = 1.0d0/makenonzero(r4m2)
1771  r5m2i = 1.0d0/makenonzero(r5m2)
1772  r6m2i = 1.0d0/makenonzero(r6m2)
1773 
1774  DO ifunnz = 1,nfunnz
1775  ifun = mapfunnz2funcorr44(ifunnz)
1776  CALL rflu_mapm2ijkl(ifun,l,m,j,k)
1777 
1778  term1 = &
1779  r1m23*r6m23*(r1m2i*(z(l ) - z(l+3)) &
1780  *(z(m ) - z(m+3)) - 4.0d0*kd(l,m)) &
1781  *(r6m2i*(z(j+6) - z(j+9)) &
1782  *(z(k+6) - z(k+9)) - 4.0d0*kd(j,k))
1783  term2 = &
1784  r2m23*r5m23*(r2m2i*(z(l ) - z(l+6)) &
1785  *(z(j ) - z(j+6)) - 4.0d0*kd(l,j)) &
1786  *(r5m2i*(z(m+3) - z(m+9)) &
1787  *(z(k+3) - z(k+9)) - 4.0d0*kd(m,k))
1788  term3 = &
1789  r3m23*r4m23*(r3m2i*(z(l ) - z(l+9)) &
1790  *(z(k ) - z(k+9)) - 4.0d0*kd(l,k)) &
1791  *(r4m2i*(z(m+3) - z(m+6)) &
1792  *(z(j+3) - z(j+6)) - 4.0d0*kd(m,j))
1793 
1794  f(ifunnz) = term1 + term2 + term3
1795  END DO ! iFunNZ
1796 
1797 ! --- End
1798 
1799  END SUBROUTINE rflu_definecorrelation542
1800 
1801 
1802 
1803 
1804 
1805 
1806 ! ******************************************************************************
1807 ! Map from vector to second-order tensor
1808 ! ******************************************************************************
1809 
1810  SUBROUTINE rflu_mapk2ij(k,i,j)
1811 
1812 ! --- Arguments
1813 
1814  INTEGER, INTENT(IN) :: k
1815  INTEGER, INTENT(OUT) :: i,j
1816 
1817 ! --- Start
1818 
1819  j = (k+2)/3 ! NOTE integer division
1820  i = k - 3*(j-1)
1821 
1822 ! --- End
1823 
1824  END SUBROUTINE rflu_mapk2ij
1825 
1826 
1827 ! ******************************************************************************
1828 ! Map from vector to third-order tensor
1829 ! ******************************************************************************
1830 
1831  SUBROUTINE rflu_mapl2ijk(l,i,j,k)
1832 
1833 ! --- Arguments
1834 
1835  INTEGER, INTENT(IN) :: l
1836  INTEGER, INTENT(OUT) :: i,j,k
1837 
1838 ! --- Start
1839 
1840  k = (l+8)/9 ! NOTE integer division
1841 
1842  j = mod((l+2)/3,3) ! NOTE integer division
1843  IF ( j == 0 ) THEN
1844  j = 3
1845  END IF ! j
1846 
1847  i = l - 3*(j-1) - 9*(k-1)
1848 
1849 ! --- End
1850 
1851  END SUBROUTINE rflu_mapl2ijk
1852 
1853 
1854 ! ******************************************************************************
1855 ! Map from vector to fourth-order tensor
1856 ! ******************************************************************************
1857 
1858  SUBROUTINE rflu_mapm2ijkl(m,i,j,k,l)
1859 
1860 ! --- Arguments
1861 
1862  INTEGER, INTENT(IN) :: m
1863  INTEGER, INTENT(OUT) :: i,j,k,l
1864 
1865 ! --- Start
1866 
1867  j = mod((m+2)/3,3) ! NOTE integer division
1868  IF ( j == 0 ) THEN
1869  j = 3
1870  END IF ! j
1871 
1872  k = mod((m+8)/9,3) ! NOTE integer division
1873  IF ( k == 0 ) THEN
1874  k = 3
1875  END IF ! k
1876 
1877  l = (m+26)/27 ! NOTE integer division
1878  i = m - 27*(l-1) - 9*(k-1) - 3*(j-1)
1879 
1880 ! --- End
1881 
1882  END SUBROUTINE rflu_mapm2ijkl
1883 
1884 
1885 
1886 ! ******************************************************************************
1887 ! Get storage position for I1
1888 ! ******************************************************************************
1889 
1890  INTEGER FUNCTION rflu_geti1posoles(l,d)
1891 
1892 ! --- Arguments
1893 
1894  INTEGER, INTENT(IN) :: d,l
1895 
1896 ! --- Start
1897 
1898  rflu_geti1posoles = l + 3*(d-1)
1899 
1900 ! --- End
1901 
1902  END FUNCTION rflu_geti1posoles
1903 
1904 
1905 ! ******************************************************************************
1906 ! Get storage position for I4
1907 ! ******************************************************************************
1908 
1909  INTEGER FUNCTION rflu_geti4posoles(l,m,d,e,nCells)
1910 
1911 ! --- Arguments
1912 
1913  INTEGER, INTENT(IN) :: d,e,l,m,ncells
1914 
1915 ! --- Start
1916 
1917  rflu_geti4posoles = m + 3*(l-1) + 9*(e-1) + 9*ncells*(d-1)
1918 
1919 ! --- End
1920 
1921  END FUNCTION rflu_geti4posoles
1922 
1923 
1924 
1925 ! ******************************************************************************
1926 ! Get storage position for L
1927 ! ******************************************************************************
1928 
1929  INTEGER FUNCTION rflu_getlposoles(j,a)
1930 
1931 ! --- Arguments
1932 
1933  INTEGER, INTENT(IN) :: a,j
1934 
1935 ! --- Start
1936 
1937  rflu_getlposoles = j + 3*(a-1)
1938 
1939 ! --- End
1940 
1941  END FUNCTION rflu_getlposoles
1942 
1943 
1944 ! ******************************************************************************
1945 ! Given storage position for L, find j and a
1946 ! ******************************************************************************
1947 
1948  SUBROUTINE rflu_getlposinvoles(loc,j,a)
1949 
1950 ! --- Arguments
1951 
1952  INTEGER, INTENT(IN) :: loc
1953  INTEGER, INTENT(OUT) :: a,j
1954 
1955 ! --- Start
1956 
1957  a = (loc+2)/3
1958  j = loc - 3*(a-1)
1959 
1960 ! --- End
1961 
1962  END SUBROUTINE rflu_getlposinvoles
1963 
1964 
1965 
1966 ! ******************************************************************************
1967 ! Get storage position for Q
1968 ! ******************************************************************************
1969 
1970  INTEGER FUNCTION rflu_getqposoles(j,k,b,g,nCells)
1971 
1972 ! --- Arguments
1973 
1974  INTEGER, INTENT(IN) :: b,g,j,k,ncells
1975 
1976 ! --- Start
1977 
1978  rflu_getqposoles = k + 3*(j-1) + 9*(g-1) + 9*ncells*(b-1)
1979 
1980 ! --- End
1981 
1982  END FUNCTION rflu_getqposoles
1983 
1984 
1985 ! ******************************************************************************
1986 ! Given storage position for Q, find j,k, and b,g
1987 ! ******************************************************************************
1988 
1989  SUBROUTINE rflu_getqposinvoles(loc,nCells,j,k,b,g)
1990 
1991 ! --- Arguments
1992 
1993  INTEGER, INTENT(IN) :: loc,ncells
1994  INTEGER, INTENT(OUT) :: b,g,j,k
1995 
1996 ! --- Start
1997 
1998  j = mod((loc+2)/3,3)
1999  IF ( j == 0 ) THEN
2000  j = 3
2001  END IF ! j
2002 
2003  g = mod((loc+8)/9,ncells)
2004  IF ( g == 0 ) THEN
2005  g = ncells
2006  END IF ! g
2007 
2008  b = (loc+9*ncells-1)/(9*ncells)
2009 
2010  k = loc - 3*(j-1) - 9*(g-1) - 9*ncells*(b-1)
2011 
2012 ! --- End
2013 
2014  END SUBROUTINE rflu_getqposinvoles
2015 
2016 
2017 
2018 
2019 
2020 ! ******************************************************************************
2021 ! Destroy geometry
2022 ! ******************************************************************************
2023 
2024  SUBROUTINE rflu_destroystencilsweightsoles(pRegion)
2025 
2026  IMPLICIT NONE
2027 
2028 ! --- parameters
2029 
2030  TYPE(t_region), POINTER :: pregion
2031 
2032 ! --- local variables
2033 
2034  INTEGER :: errorflag
2035  TYPE(t_grid), POINTER :: pgrid
2036  TYPE(t_global), POINTER :: global
2037 
2038 ! ==============================================================================
2039 ! Start
2040 ! ==============================================================================
2041 
2042  global => pregion%global
2043 
2044  CALL registerfunction(global,'RFLU_DestroyStencilsWeightsOLES',&
2045  'RFLU_ModOLES.F90')
2046 
2047  IF ( global%myProcid == masterproc .AND. &
2048  global%verbLevel > verbose_none ) THEN
2049  WRITE(stdout,'(A)') solver_name
2050  WRITE(stdout,'(A,1X,A)') solver_name,'Destroying optimal LES '// &
2051  'stencils and weights...'
2052  END IF ! global%verbLevel
2053 
2054 ! ==============================================================================
2055 ! Set grid pointer
2056 ! ==============================================================================
2057 
2058  pgrid => pregion%grid
2059 
2060 ! ==============================================================================
2061 ! Deallocate memory
2062 ! ==============================================================================
2063 
2064  DEALLOCATE(pgrid%fsOLES,stat=errorflag)
2065  global%error = errorflag
2066  IF ( global%error /= err_none ) THEN
2067  CALL errorstop(global,err_deallocate,__line__,'region%grid%fsOLES')
2068  END IF ! global%error
2069 
2070  DEALLOCATE(pgrid%intLimOLES,stat=errorflag)
2071  global%error = errorflag
2072  IF ( global%error /= err_none ) THEN
2073  CALL errorstop(global,err_deallocate,__line__,'region%grid%intLimOLES')
2074  END IF ! global%error
2075 
2076 ! ==============================================================================
2077 ! End
2078 ! ==============================================================================
2079 
2080  IF ( global%myProcid == masterproc .AND. &
2081  global%verbLevel > verbose_none ) THEN
2082  WRITE(stdout,'(A,1X,A)') solver_name,'Destroying optimal LES '// &
2083  'stencils and weights done.'
2084  END IF ! global%verbLevel
2085 
2086  CALL deregisterfunction(global)
2087 
2088  END SUBROUTINE rflu_destroystencilsweightsoles
2089 
2090 
2091 
2092 
2093 ! ******************************************************************************
2094 ! Allocate DCUHRE arrays
2095 ! ******************************************************************************
2096 
2097  SUBROUTINE rflu_allocatedcuhrearrays(global,nDim,nFunNZ,nFun)
2098 
2099  IMPLICIT NONE
2100 
2101 ! --- parameters
2102 
2103  INTEGER, INTENT(IN) :: ndim,nfunnz,nfun
2104  TYPE(t_global), POINTER :: global
2105 
2106 ! --- locals
2107 
2108  INTEGER :: errorflag
2109 
2110 ! ==============================================================================
2111 ! Start
2112 ! ==============================================================================
2113 
2114  CALL registerfunction(global,'RFLU_AllocateDCUHREArrays',&
2115  'RFLU_ModOLES.F90')
2116 
2117 ! ==============================================================================
2118 ! Allocate memory
2119 ! ==============================================================================
2120 
2121  ALLOCATE(lowlim(ndim),stat=errorflag)
2122  global%error = errorflag
2123  IF ( global%error /= err_none ) THEN
2124  CALL errorstop(global,err_allocate,__line__,'lowLim')
2125  END IF ! global%error
2126 
2127  ALLOCATE(upplim(ndim),stat=errorflag)
2128  global%error = errorflag
2129  IF ( global%error /= err_none ) THEN
2130  CALL errorstop(global,err_allocate,__line__,'uppLim')
2131  END IF ! global%error
2132 
2133  ALLOCATE(errabsest(nfun),stat=errorflag)
2134  global%error = errorflag
2135  IF ( global%error /= err_none ) THEN
2136  CALL errorstop(global,err_allocate,__line__,'errAbsEst')
2137  END IF ! global%error
2138 
2139  ALLOCATE(integralnz(nfunnz),stat=errorflag)
2140  global%error = errorflag
2141  IF ( global%error /= err_none ) THEN
2142  CALL errorstop(global,err_allocate,__line__,'integralNZ')
2143  END IF ! global%error
2144 
2145 ! ==============================================================================
2146 ! End
2147 ! ==============================================================================
2148 
2149  CALL deregisterfunction(global)
2150 
2151  END SUBROUTINE rflu_allocatedcuhrearrays
2152 
2153 
2154 
2155 
2156 
2157 ! ******************************************************************************
2158 ! Deallocate DCUHRE arrays
2159 ! ******************************************************************************
2160 
2161  SUBROUTINE rflu_deallocatedcuhrearrays(global)
2162 
2163  IMPLICIT NONE
2164 
2165 ! --- arguments
2166 
2167  TYPE(t_global), POINTER :: global
2168 
2169 ! --- locals
2170 
2171  INTEGER :: errorflag
2172 
2173 ! ==============================================================================
2174 ! Start
2175 ! ==============================================================================
2176 
2177  CALL registerfunction(global,'RFLU_DeallocateDCUHREArrays',&
2178  'RFLU_ModOLES.F90')
2179 
2180 ! ==============================================================================
2181 ! Allocate memory
2182 ! ==============================================================================
2183 
2184  DEALLOCATE(lowlim,stat=errorflag)
2185  global%error = errorflag
2186  IF ( global%error /= err_none ) THEN
2187  CALL errorstop(global,err_deallocate,__line__,'lowLim')
2188  END IF ! global%error
2189 
2190  DEALLOCATE(upplim,stat=errorflag)
2191  global%error = errorflag
2192  IF ( global%error /= err_none ) THEN
2193  CALL errorstop(global,err_deallocate,__line__,'uppLim')
2194  END IF ! global%error
2195 
2196  DEALLOCATE(errabsest,stat=errorflag)
2197  global%error = errorflag
2198  IF ( global%error /= err_none ) THEN
2199  CALL errorstop(global,err_deallocate,__line__,'errAbsEst')
2200  END IF ! global%error
2201 
2202  DEALLOCATE(integralnz,stat=errorflag)
2203  global%error = errorflag
2204  IF ( global%error /= err_none ) THEN
2205  CALL errorstop(global,err_deallocate,__line__,'integralNZ')
2206  END IF ! global%error
2207 
2208 ! ==============================================================================
2209 ! End
2210 ! ==============================================================================
2211 
2212  CALL deregisterfunction(global)
2213 
2214  END SUBROUTINE rflu_deallocatedcuhrearrays
2215 
2216 
2217 
2218 
2219 ! ******************************************************************************
2220 ! Set mapping for second-order two-point correlation (for I^2)
2221 ! ******************************************************************************
2222 
2223  SUBROUTINE rflu_setmapfunnz2funcorr22(nFunNZ)
2224 
2225 ! --- Arguments
2226 
2227  INTEGER, INTENT(IN) :: nfunnz
2228 
2229 ! --- Locals
2230 
2231  INTEGER :: ifun,ifunnz,j,l
2232  INTEGER, PARAMETER :: nfun = 9
2233 
2234 ! --- Start
2235 
2236  ifunnz = 0
2237 
2238  IF ( nfunnz /= nfun ) THEN
2239  DO ifun = 1,nfun
2240  CALL rflu_mapk2ij(ifun,l,j)
2241 
2242  IF ( l == j ) THEN
2243  ifunnz = ifunnz + 1
2244  mapfunnz2funcorr22(ifunnz) = ifun
2245  END IF ! i
2246  END DO ! iFun
2247 
2248  mapfunnz2funcorr22(ifunnz+1:nfun) = crazy_value_int
2249  ELSE
2250  DO ifun = 1,nfun
2251  ifunnz = ifunnz + 1
2252  mapfunnz2funcorr22(ifunnz) = ifun
2253  END DO ! iFun
2254  END IF ! nFunNZ
2255 
2256 ! --- End
2257 
2258  END SUBROUTINE rflu_setmapfunnz2funcorr22
2259 
2260 
2261 
2262 ! ******************************************************************************
2263 ! Set mapping for third-order two-point correlation (for I^1)
2264 ! ******************************************************************************
2265 
2266  SUBROUTINE rflu_setmapfunnz2funcorr32(nFunNZ)
2267 
2268 ! --- Arguments
2269 
2270  INTEGER, INTENT(IN) :: nfunnz
2271 
2272 ! --- Locals
2273 
2274  INTEGER :: i,ifun,ifunnz,l
2275  INTEGER, PARAMETER :: nfun = 9
2276 
2277 ! --- Start
2278 
2279  ifunnz = 0
2280 
2281  IF ( nfunnz /= nfun ) THEN
2282  DO ifun = 1,nfun
2283  CALL rflu_mapk2ij(ifun,l,i)
2284 
2285  IF ( l == i ) THEN
2286  ifunnz = ifunnz + 1
2287  mapfunnz2funcorr32(ifunnz) = ifun
2288  END IF ! i
2289  END DO ! iFun
2290 
2291  mapfunnz2funcorr32(ifunnz+1:nfun) = crazy_value_int
2292  ELSE
2293  DO ifun = 1,nfun
2294  ifunnz = ifunnz + 1
2295  mapfunnz2funcorr32(ifunnz) = ifun
2296  END DO ! iFun
2297  END IF ! nFunNZ
2298 
2299 ! --- End
2300 
2301  END SUBROUTINE rflu_setmapfunnz2funcorr32
2302 
2303 
2304 
2305 
2306 ! ******************************************************************************
2307 ! Set mapping for fourth-order three-point correlation (for I^4)
2308 ! ******************************************************************************
2309 
2310  SUBROUTINE rflu_setmapfunnz2funcorr43(nFunNZ)
2311 
2312 ! --- Arguments
2313 
2314  INTEGER, INTENT(IN) :: nfunnz
2315 
2316 ! --- Locals
2317 
2318  INTEGER :: i,ifun,ifunnz,l,m,s
2319  INTEGER, PARAMETER :: nfun = 27
2320  REAL(RFREAL) :: term
2321 
2322 ! --- Start
2323 
2324  ifunnz = 0
2325 
2326  s = nzloc
2327 
2328  IF ( nfunnz /= nfun ) THEN
2329  DO ifun = 1,nfun
2330  CALL rflu_mapl2ijk(ifun,l,m,i)
2331 
2332  term = kd(l,m)*kd(i,s) + kd(l,i)*kd(m,s) + kd(l,s)*kd(m,i)
2333 
2334  IF ( nint(term) /= 0 ) THEN
2335  ifunnz = ifunnz + 1
2336  mapfunnz2funcorr43(ifunnz) = ifun
2337  END IF ! i
2338  END DO ! iFun
2339 
2340  mapfunnz2funcorr43(ifunnz+1:nfun) = crazy_value_int
2341  ELSE
2342  DO ifun = 1,nfun
2343  ifunnz = ifunnz + 1
2344  mapfunnz2funcorr43(ifunnz) = ifun
2345  END DO ! iFun
2346  END IF ! nFunNZ
2347 
2348 ! --- End
2349 
2350  END SUBROUTINE rflu_setmapfunnz2funcorr43
2351 
2352 
2353 
2354 
2355 ! ******************************************************************************
2356 ! Set mapping for fourth-order four-point correlation (for I^5)
2357 ! ******************************************************************************
2358 
2359  SUBROUTINE rflu_setmapfunnz2funcorr44(nFunNZ)
2360 
2361 ! --- Arguments
2362 
2363  INTEGER, INTENT(IN) :: nfunnz
2364 
2365 ! --- Locals
2366 
2367  INTEGER :: ifun,ifunnz,j,k,l,m
2368  INTEGER, PARAMETER :: nfun = 81
2369  REAL(RFREAL) :: term
2370 
2371 ! --- Start
2372 
2373  ifunnz = 0
2374 
2375  IF ( nfunnz /= nfun ) THEN
2376  DO ifun = 1,nfun
2377  CALL rflu_mapm2ijkl(ifun,l,m,j,k)
2378 
2379  term = kd(l,m)*kd(j,k) + kd(l,j)*kd(m,k) + kd(l,k)*kd(m,j)
2380 
2381  IF ( nint(term) /= 0 ) THEN
2382  ifunnz = ifunnz + 1
2383  mapfunnz2funcorr44(ifunnz) = ifun
2384  END IF ! i
2385  END DO ! iFun
2386 
2387  mapfunnz2funcorr44(ifunnz+1:nfun) = crazy_value_int
2388  ELSE
2389  DO ifun = 1,nfun
2390  ifunnz = ifunnz + 1
2391  mapfunnz2funcorr44(ifunnz) = ifun
2392  END DO ! iFun
2393  END IF ! nFunNZ
2394 
2395 ! --- End
2396 
2397  END SUBROUTINE rflu_setmapfunnz2funcorr44
2398 
2399 
2400 
2401 
2402 ! ******************************************************************************
2403 ! End
2404 ! ******************************************************************************
2405 
2406 END MODULE rflu_modoles
2407 
2408 
2409 ! ******************************************************************************
2410 !
2411 ! RCS Revision history:
2412 !
2413 ! $Log: RFLU_ModOLES.F90,v $
2414 ! Revision 1.10 2008/12/06 08:44:23 mtcampbe
2415 ! Updated license.
2416 !
2417 ! Revision 1.9 2008/11/19 22:17:34 mtcampbe
2418 ! Added Illinois Open Source License/Copyright
2419 !
2420 ! Revision 1.8 2004/01/22 16:03:59 haselbac
2421 ! Made contents of modules PRIVATE, only procs PUBLIC, to avoid errors on ALC and titan
2422 !
2423 ! Revision 1.7 2004/01/11 02:17:34 jiao
2424 ! Eliminated some redundant trailing spaces that made some lines too long.
2425 ! This changed was necessary to compile with NAG F90 compiler.
2426 !
2427 ! Revision 1.6 2003/05/16 02:27:44 haselbac
2428 ! Removed KIND=RFREAL from NINT statements
2429 !
2430 ! Revision 1.5 2003/03/15 18:14:17 haselbac
2431 ! Added KIND qualifyer to NINT statements
2432 !
2433 ! Revision 1.4 2002/10/08 15:49:21 haselbac
2434 ! {IO}STAT=global%error replaced by {IO}STAT=errorFlag - SGI problem
2435 !
2436 ! Revision 1.3 2002/09/09 15:10:35 haselbac
2437 ! global now under regions, many bug fixes as a result of CTR stay
2438 !
2439 ! Revision 1.2 2002/07/27 18:11:28 haselbac
2440 ! Corrected I41, I42, I51, and I52 correlations, now get correct symmetries
2441 !
2442 ! Revision 1.1 2002/07/25 14:51:54 haselbac
2443 ! Initial revision
2444 !
2445 ! ******************************************************************************
2446 
2447 
2448 
2449 
2450 
2451 
2452 
2453 
2454 
2455 
2456 
2457 
2458 
2459 
2460 
2461 
Aff_transformation_2::m
FT m(int i, int j) const
Definition: Aff_transformation_2.h:228
rflu_modoles::rflu_definecorrelation221
subroutine, public rflu_definecorrelation221(nDim, z, nFunNZ, f)
Definition: RFLU_ModOLES.F90:1333
rflu_modoles::rflu_computegeometrictermsoles
subroutine, public rflu_computegeometrictermsoles(pRegion)
Definition: RFLU_ModOLES.F90:881
Bbox_2::ymin
double ymin() const
Definition: CGAL/include/CGAL/Bbox_2.h:101
d
const NT & d
Definition: rational_rotation.h:73
Bbox_2::xmax
double xmax() const
Definition: CGAL/include/CGAL/Bbox_2.h:106
rflu_modoctree::rflu_queryoctree
subroutine, public rflu_queryoctree(XPT, YPT, ZPT, NUMP, NEIGHP)
Definition: RFLU_ModOctree.F90:405
k
j indices k indices k
Definition: Indexing.h:6
rflu_modoles::rflu_definecorrelation432
subroutine, public rflu_definecorrelation432(nDim, z, nFunNZ, f)
Definition: RFLU_ModOLES.F90:1536
s
double s
Definition: blastest.C:80
rflu_modoles::rflu_definecorrelation430
subroutine, public rflu_definecorrelation430(nDim, z, nFunNZ, f)
Definition: RFLU_ModOLES.F90:1415
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
Bbox_2::xmin
double xmin() const
Definition: CGAL/include/CGAL/Bbox_2.h:96
moderror::registerfunction
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
Color::b
unsigned char b() const
Definition: Color.h:70
ic
**********************************************************************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 ic
Definition: RFLO_FindSourceCell.h:27
Bbox_3::zmin
double zmin() const
Definition: CGAL/include/CGAL/Bbox_3.h:106
modgrid::t_grid
Definition: ModGrid.F90:56
sqrt
double sqrt(double d)
Definition: double.h:73
rflu_modoctree::rflu_buildoctree
subroutine, public rflu_buildoctree(XI, YI, ZI, XLOW, XUPP, YLOW, YUPP, ZLOW, ZUPP)
Definition: RFLU_ModOctree.F90:152
rflu_modoles::rflu_geti1posoles
INTEGER function, public rflu_geti1posoles(l, d)
Definition: RFLU_ModOLES.F90:1890
rflu_modoles::rflu_getqposinvoles
subroutine, public rflu_getqposinvoles(loc, nCells, j, k, b, g)
Definition: RFLU_ModOLES.F90:1989
rflu_modoles::rflu_setmapfunnz2funcorr32
subroutine, public rflu_setmapfunnz2funcorr32(nFunNZ)
Definition: RFLU_ModOLES.F90:2266
rflu_allocatedcuhrearrays
subroutine rflu_allocatedcuhrearrays(nDim, nFunNZ, nFun)
Definition: RFLU_AllocateDCUHREArrays.F90:45
rflu_modoles::rflu_findprototypefacesoles
subroutine, public rflu_findprototypefacesoles(pRegion)
Definition: RFLU_ModOLES.F90:742
rflu_modoles::rflu_setmapfunnz2funcorr44
subroutine, public rflu_setmapfunnz2funcorr44(nFunNZ)
Definition: RFLU_ModOLES.F90:2359
rflu_modoles
Definition: RFLU_ModOLES.F90:39
rflu_modoles::rflu_getlposoles
INTEGER function, public rflu_getlposoles(j, a)
Definition: RFLU_ModOLES.F90:1929
rflu_modoles::rflu_getlposinvoles
subroutine, public rflu_getlposinvoles(loc, j, a)
Definition: RFLU_ModOLES.F90:1948
rflu_modoles::rflu_definecorrelation542
subroutine, public rflu_definecorrelation542(nDim, z, nFunNZ, f)
Definition: RFLU_ModOLES.F90:1724
rflu_modoles::rflu_definecorrelation431
subroutine, public rflu_definecorrelation431(nDim, z, nFunNZ, f)
Definition: RFLU_ModOLES.F90:1450
rflu_modoles::rflu_enforcesymmetryoles
subroutine, public rflu_enforcesymmetryoles(pRegion)
Definition: RFLU_ModOLES.F90:1179
z
void int int int REAL REAL REAL * z
Definition: write.cpp:76
rflu_modoles::rflu_mapk2ij
subroutine, public rflu_mapk2ij(k, i, j)
Definition: RFLU_ModOLES.F90:1810
rflu_modoctree::rflu_destroyoctree
subroutine, public rflu_destroyoctree(global)
Definition: RFLU_ModOctree.F90:668
Mesh::GenericElement::nfaces
IndexType nfaces() const
Definition: Mesh.H:641
rflu_modoles::rflu_definecorrelation32
subroutine, public rflu_definecorrelation32(nDim, z, nFunNZ, f)
Definition: RFLU_ModOLES.F90:1369
rflu_modoles::rflu_createintegralsoles
subroutine, public rflu_createintegralsoles(pRegion)
Definition: RFLU_ModOLES.F90:283
rflu_modoles::rflu_definecorrelation540
subroutine, public rflu_definecorrelation540(nDim, z, nFunNZ, f)
Definition: RFLU_ModOLES.F90:1612
i
blockLoc i
Definition: read.cpp:79
Bbox_3::zmax
double zmax() const
Definition: CGAL/include/CGAL/Bbox_3.h:121
rflu_modoles::rflu_setmapfunnz2funcorr22
subroutine, public rflu_setmapfunnz2funcorr22(nFunNZ)
Definition: RFLU_ModOLES.F90:2223
n
const NT & n
Definition: rational_rotation.h:72
rflu_modoles::rflu_mapm2ijkl
subroutine, public rflu_mapm2ijkl(m, i, j, k, l)
Definition: RFLU_ModOLES.F90:1858
Bbox_2::ymax
double ymax() const
Definition: CGAL/include/CGAL/Bbox_2.h:111
rflu_deallocatedcuhrearrays
subroutine rflu_deallocatedcuhrearrays
Definition: RFLU_AllocateDCUHREArrays.F90:108
rflu_modoles::rflu_definecorrelation541
subroutine, public rflu_definecorrelation541(nDim, z, nFunNZ, f)
Definition: RFLU_ModOLES.F90:1644
modsortsearch::quicksortrfrealinteger
subroutine quicksortrfrealinteger(a, b, n)
Definition: ModSortSearch.F90:840
min
Vector_n min(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:346
rflu_modoles::rflu_geti4posoles
INTEGER function, public rflu_geti4posoles(l, m, d, e, nCells)
Definition: RFLU_ModOLES.F90:1909
rflu_modoctree
Definition: RFLU_ModOctree.F90:51
j
j indices j
Definition: Indexing.h:6
rflu_modoles::rflu_definecorrelation220
subroutine, public rflu_definecorrelation220(nDim, z, nFunNZ, f)
Definition: RFLU_ModOLES.F90:1304
rflu_modoles::rflu_buildsymmetrymapsoles
subroutine, public rflu_buildsymmetrymapsoles(pRegion)
Definition: RFLU_ModOLES.F90:1106
modglobal::t_global
Definition: ModGlobal.F90:50
moderror::errorstop
subroutine errorstop(global, errorCode, errorLine, addMessage)
Definition: ModError.F90:483
rflu_modoles::rflu_buildstencilsoles
subroutine, public rflu_buildstencilsoles(pRegion)
Definition: RFLU_ModOLES.F90:483
dist
long double dist(long double *coord1, long double *coord2, int size)
Definition: utilities/rc_cmp/util.cpp:305
rflu_modoles::rflu_mapl2ijk
subroutine, public rflu_mapl2ijk(l, i, j, k)
Definition: RFLU_ModOLES.F90:1831
rflu_modoles::rflu_getqposoles
INTEGER function, public rflu_getqposoles(j, k, b, g, nCells)
Definition: RFLU_ModOLES.F90:1970
moderror::deregisterfunction
subroutine deregisterfunction(global)
Definition: ModError.F90:469
rflu_modoctree::rflu_createoctree
subroutine, public rflu_createoctree(global, nPoints)
Definition: RFLU_ModOctree.F90:92
rflu_modoles::rflu_destroystencilsweightsoles
subroutine, public rflu_destroystencilsweightsoles(pRegion)
Definition: RFLU_ModOLES.F90:2024
rflu_modoles::rflu_createstencilsweightsoles
subroutine, public rflu_createstencilsweightsoles(pRegion)
Definition: RFLU_ModOLES.F90:149
modtools::floatequal
LOGICAL function floatequal(a, b, tolIn)
Definition: ModTools.F90:99
rflu_modoles::rflu_setmapfunnz2funcorr43
subroutine, public rflu_setmapfunnz2funcorr43(nFunNZ)
Definition: RFLU_ModOLES.F90:2310
Line_2::a
RT a() const
Definition: Line_2.h:140
Color::g
unsigned char g() const
Definition: Color.h:69
modtools::makenonzero
real(rfreal) function makenonzero(x)
Definition: ModTools.F90:85
modbndpatch::t_patch
Definition: ModBndPatch.F90:91