Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLU_ComputeGridSpeeds.F90
Go to the documentation of this file.
1 ! *********************************************************************
2 ! * Rocstar Simulation Suite *
3 ! * Copyright@2015, Illinois Rocstar LLC. All rights reserved. *
4 ! * *
5 ! * Illinois Rocstar LLC *
6 ! * Champaign, IL *
7 ! * www.illinoisrocstar.com *
8 ! * sales@illinoisrocstar.com *
9 ! * *
10 ! * License: See LICENSE file in top level of distribution package or *
11 ! * http://opensource.org/licenses/NCSA *
12 ! *********************************************************************
13 ! *********************************************************************
14 ! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
15 ! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
16 ! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
17 ! * NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
18 ! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
19 ! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
20 ! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
21 ! * USE OR OTHER DEALINGS WITH THE SOFTWARE. *
22 ! *********************************************************************
23 ! ******************************************************************************
24 !
25 ! Purpose: Compute grid speeds.
26 !
27 ! Description: The grid speeds are computed by computing the volume swept out
28 ! by the motion of the face, divided by a surface area and the time step.
29 ! Depending on whether the face is triangular or quadrilateral, the volume
30 ! swept out by the face motion is a prism or a hexahedron.
31 !
32 ! Input:
33 ! pRegion Pointer to region
34 !
35 ! Output: None.
36 !
37 ! Notes:
38 ! 1. The order in which RFLU_ComputeGridSpeeds and RFLU_BuildGeometry are
39 ! called is important because the grid speeds are defined to be a volume
40 ! change divided by a face area (and the time-step). If the face area is
41 ! from the old grid, the grid speeds will not satisfy continuity when
42 ! used to compute fluxes.
43 ! 2. The volumes of the polyhedra swept out by the triangular and
44 ! quadrilateral faces are computed in a manner consistent with the way in
45 ! which the cell volumes are computed, see, e.g., RFLU_ModGeometry.F90.
46 !
47 ! *****************************************************************************
48 !
49 ! $Id: RFLU_ComputeGridSpeeds.F90,v 1.9 2008/12/06 08:44:29 mtcampbe Exp $
50 !
51 ! Copyright: (c) 2002-2005 by the University of Illinois
52 !
53 ! ******************************************************************************
54 
55 SUBROUTINE rflu_computegridspeeds(pRegion)
56 
57  USE moddatatypes
58  USE modgrid, ONLY: t_grid
59  USE moddatastruct, ONLY: t_region
60  USE modbndpatch, ONLY: t_patch
61  USE modglobal, ONLY: t_global
62  USE moderror
63  USE modmpi
64  USE modparameters
65 
66  USE rflu_modgrid
67 
68  IMPLICIT NONE
69 
70 ! ******************************************************************************
71 ! Declarations and definitions
72 ! ******************************************************************************
73 
74 ! ==============================================================================
75 ! Arguments
76 ! ==============================================================================
77 
78  TYPE(t_region), POINTER :: pregion
79 
80 ! ==============================================================================
81 ! Locals
82 ! ==============================================================================
83 
84  CHARACTER(CHRLEN) :: rcsidentstring
85  INTEGER :: errorflag,ifc,ifcl,ipatch,iv,v1,v1l,v2,v2l,v3,v3l,v4,v4l
86  REAL(RFREAL), PARAMETER :: thrd = 1.0_rfreal/3.0_rfreal
87  REAL(RFREAL) :: fcx,fcy,fcz,fnx,fny,fnz,term
88  REAL(RFREAL) :: xyzavg(xcoord:zcoord)
89  REAL(RFREAL) :: xyzhex(xcoord:zcoord,8),xyznodes(xcoord:zcoord,4), &
90  xyzpri(xcoord:zcoord,6)
91  REAL(RFREAL) :: vsnew, vsold
92  REAL(RFREAL), DIMENSION(:,:), POINTER :: pxyz,pxyzold,pxyzold2
93  TYPE(t_global), POINTER :: global
94  TYPE(t_grid), POINTER :: pgrid,pgridold,pgridold2
95  TYPE(t_patch), POINTER :: ppatch
96 
97 ! ******************************************************************************
98 ! Start
99 ! ******************************************************************************
100 
101  rcsidentstring = '$RCSfile: RFLU_ComputeGridSpeeds.F90,v $ $Revision: 1.9 $'
102 
103  global => pregion%global
104 
105  CALL registerfunction(global,'RFLU_ComputeGridSpeeds',&
106  'RFLU_ComputeGridSpeeds.F90')
107 
108  IF ( global%myProcid == masterproc .AND. &
109  global%verbLevel >= verbose_high ) THEN
110  WRITE(stdout,'(A,1X,A)') solver_name,'Computing grid speeds...'
111  WRITE(stdout,'(A,3X,A)') solver_name,'Grid speed extrema:'
112  END IF ! global%myProcid
113 
114 ! ******************************************************************************
115 ! Set pointers and variables
116 ! ******************************************************************************
117 
118  pgrid => pregion%grid
119  pgridold => pregion%gridOld
120  pgridold2 => pregion%gridOld2
121 
122  pxyz => pgrid%xyz
123  pxyzold => pgridold%xyz
124  pxyzold2 => pgridold2%xyz
125 
126 ! ******************************************************************************
127 ! Interior faces
128 ! ******************************************************************************
129 
130  DO ifc = 1,pgrid%nFaces
131  v1 = pgrid%f2v(1,ifc)
132  v2 = pgrid%f2v(2,ifc)
133  v3 = pgrid%f2v(3,ifc)
134  v4 = pgrid%f2v(4,ifc)
135 
136  pgrid%gs(ifc) = 0.0_rfreal
137 
138  IF ( global%solverType == solv_implicit_nk ) THEN
139  vsnew = 0.0_rfreal
140  vsold = 0.0_rfreal
141  END IF ! global%solverType
142 
143 ! ==============================================================================
144 ! Triangular face, sweeps out prism during grid motion
145 ! ==============================================================================
146 
147  IF ( v4 == vert_none ) THEN ! triangular face
148  xyzpri(xcoord,1) = pxyzold(xcoord,v1)
149  xyzpri(xcoord,2) = pxyzold(xcoord,v2)
150  xyzpri(xcoord,3) = pxyzold(xcoord,v3)
151 
152  xyzpri(ycoord,1) = pxyzold(ycoord,v1)
153  xyzpri(ycoord,2) = pxyzold(ycoord,v2)
154  xyzpri(ycoord,3) = pxyzold(ycoord,v3)
155 
156  xyzpri(zcoord,1) = pxyzold(zcoord,v1)
157  xyzpri(zcoord,2) = pxyzold(zcoord,v2)
158  xyzpri(zcoord,3) = pxyzold(zcoord,v3)
159 
160  xyzpri(xcoord,4) = pxyz(xcoord,v1)
161  xyzpri(xcoord,5) = pxyz(xcoord,v2)
162  xyzpri(xcoord,6) = pxyz(xcoord,v3)
163 
164  xyzpri(ycoord,4) = pxyz(ycoord,v1)
165  xyzpri(ycoord,5) = pxyz(ycoord,v2)
166  xyzpri(ycoord,6) = pxyz(ycoord,v3)
167 
168  xyzpri(zcoord,4) = pxyz(zcoord,v1)
169  xyzpri(zcoord,5) = pxyz(zcoord,v2)
170  xyzpri(zcoord,6) = pxyz(zcoord,v3)
171 
172 ! ------------------------------------------------------------------------------
173 ! Compute average coordinates (approximate centroid)
174 ! ------------------------------------------------------------------------------
175 
176  xyzavg(xcoord) = 0.0_rfreal
177  xyzavg(ycoord) = 0.0_rfreal
178  xyzavg(zcoord) = 0.0_rfreal
179 
180  term = 1.0_rfreal/6.0_rfreal
181 
182  DO iv = 1,6
183  xyzavg(xcoord) = xyzavg(xcoord) + xyzpri(xcoord,iv)
184  xyzavg(ycoord) = xyzavg(ycoord) + xyzpri(ycoord,iv)
185  xyzavg(zcoord) = xyzavg(zcoord) + xyzpri(zcoord,iv)
186  END DO ! iv
187 
188  xyzavg(xcoord) = term*xyzavg(xcoord)
189  xyzavg(ycoord) = term*xyzavg(ycoord)
190  xyzavg(zcoord) = term*xyzavg(zcoord)
191 
192 ! ------------------------------------------------------------------------------
193 ! Loop over the five faces of the prism to compute volume
194 ! ------------------------------------------------------------------------------
195 
196  DO ifcl = 1,5
197  v1l = f2vpri(1,ifcl)
198  v2l = f2vpri(2,ifcl)
199  v3l = f2vpri(3,ifcl)
200  v4l = f2vpri(4,ifcl)
201 
202  IF ( v4l == 0 ) THEN ! triangular face
203  xyznodes(xcoord,1) = xyzpri(xcoord,v1l) - xyzavg(xcoord)
204  xyznodes(xcoord,2) = xyzpri(xcoord,v2l) - xyzavg(xcoord)
205  xyznodes(xcoord,3) = xyzpri(xcoord,v3l) - xyzavg(xcoord)
206 
207  xyznodes(ycoord,1) = xyzpri(ycoord,v1l) - xyzavg(ycoord)
208  xyznodes(ycoord,2) = xyzpri(ycoord,v2l) - xyzavg(ycoord)
209  xyznodes(ycoord,3) = xyzpri(ycoord,v3l) - xyzavg(ycoord)
210 
211  xyznodes(zcoord,1) = xyzpri(zcoord,v1l) - xyzavg(zcoord)
212  xyznodes(zcoord,2) = xyzpri(zcoord,v2l) - xyzavg(zcoord)
213  xyznodes(zcoord,3) = xyzpri(zcoord,v3l) - xyzavg(zcoord)
214 
215  CALL facevectortria(xyznodes(xcoord:zcoord,1:3),fnx,fny,fnz)
216  CALL facecentroidtria(xyznodes(xcoord:zcoord,1:3),fcx,fcy,fcz)
217  ELSE
218  xyznodes(xcoord,1) = xyzpri(xcoord,v1l) - xyzavg(xcoord)
219  xyznodes(xcoord,2) = xyzpri(xcoord,v2l) - xyzavg(xcoord)
220  xyznodes(xcoord,3) = xyzpri(xcoord,v3l) - xyzavg(xcoord)
221  xyznodes(xcoord,4) = xyzpri(xcoord,v4l) - xyzavg(xcoord)
222 
223  xyznodes(ycoord,1) = xyzpri(ycoord,v1l) - xyzavg(ycoord)
224  xyznodes(ycoord,2) = xyzpri(ycoord,v2l) - xyzavg(ycoord)
225  xyznodes(ycoord,3) = xyzpri(ycoord,v3l) - xyzavg(ycoord)
226  xyznodes(ycoord,4) = xyzpri(ycoord,v4l) - xyzavg(ycoord)
227 
228  xyznodes(zcoord,1) = xyzpri(zcoord,v1l) - xyzavg(zcoord)
229  xyznodes(zcoord,2) = xyzpri(zcoord,v2l) - xyzavg(zcoord)
230  xyznodes(zcoord,3) = xyzpri(zcoord,v3l) - xyzavg(zcoord)
231  xyznodes(zcoord,4) = xyzpri(zcoord,v4l) - xyzavg(zcoord)
232 
233  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4),fnx,fny,fnz)
234  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4),fcx,fcy,fcz)
235  END IF ! v4l
236 
237  IF ( global%solverType == solv_implicit_nk ) THEN
238  vsnew = vsnew + (fcx*fnx + fcy*fny + fcz*fnz)
239  ELSE
240  pgrid%gs(ifc) = pgrid%gs(ifc) + (fcx*fnx + fcy*fny + fcz*fnz)
241  END IF ! global%solverType
242 
243  END DO ! ifcl
244 
245 ! ==============================================================================
246 ! Quadrilateral face, sweeps out hexahedron during grid motion
247 ! ==============================================================================
248 
249  ELSE ! quadrilateral face
250  xyzhex(xcoord,1) = pxyzold(xcoord,v1)
251  xyzhex(xcoord,2) = pxyzold(xcoord,v2)
252  xyzhex(xcoord,3) = pxyzold(xcoord,v3)
253  xyzhex(xcoord,4) = pxyzold(xcoord,v4)
254 
255  xyzhex(ycoord,1) = pxyzold(ycoord,v1)
256  xyzhex(ycoord,2) = pxyzold(ycoord,v2)
257  xyzhex(ycoord,3) = pxyzold(ycoord,v3)
258  xyzhex(ycoord,4) = pxyzold(ycoord,v4)
259 
260  xyzhex(zcoord,1) = pxyzold(zcoord,v1)
261  xyzhex(zcoord,2) = pxyzold(zcoord,v2)
262  xyzhex(zcoord,3) = pxyzold(zcoord,v3)
263  xyzhex(zcoord,4) = pxyzold(zcoord,v4)
264 
265  xyzhex(xcoord,5) = pxyz(xcoord,v1)
266  xyzhex(xcoord,6) = pxyz(xcoord,v2)
267  xyzhex(xcoord,7) = pxyz(xcoord,v3)
268  xyzhex(xcoord,8) = pxyz(xcoord,v4)
269 
270  xyzhex(ycoord,5) = pxyz(ycoord,v1)
271  xyzhex(ycoord,6) = pxyz(ycoord,v2)
272  xyzhex(ycoord,7) = pxyz(ycoord,v3)
273  xyzhex(ycoord,8) = pxyz(ycoord,v4)
274 
275  xyzhex(zcoord,5) = pxyz(zcoord,v1)
276  xyzhex(zcoord,6) = pxyz(zcoord,v2)
277  xyzhex(zcoord,7) = pxyz(zcoord,v3)
278  xyzhex(zcoord,8) = pxyz(zcoord,v4)
279 
280 ! ------------------------------------------------------------------------------
281 ! Compute average coordinates (approximate centroid)
282 ! ------------------------------------------------------------------------------
283 
284  xyzavg(xcoord) = 0.0_rfreal
285  xyzavg(ycoord) = 0.0_rfreal
286  xyzavg(zcoord) = 0.0_rfreal
287 
288  term = 1.0_rfreal/8.0_rfreal
289 
290  DO iv = 1,8
291  xyzavg(xcoord) = xyzavg(xcoord) + xyzhex(xcoord,iv)
292  xyzavg(ycoord) = xyzavg(ycoord) + xyzhex(ycoord,iv)
293  xyzavg(zcoord) = xyzavg(zcoord) + xyzhex(zcoord,iv)
294  END DO ! iv
295 
296  xyzavg(xcoord) = term*xyzavg(xcoord)
297  xyzavg(ycoord) = term*xyzavg(ycoord)
298  xyzavg(zcoord) = term*xyzavg(zcoord)
299 
300 ! ------------------------------------------------------------------------------
301 ! Loop over the six faces of the prism to compute volume
302 ! ------------------------------------------------------------------------------
303 
304  DO ifcl = 1,6
305  v1l = f2vhex(1,ifcl)
306  v2l = f2vhex(2,ifcl)
307  v3l = f2vhex(3,ifcl)
308  v4l = f2vhex(4,ifcl)
309 
310  xyznodes(xcoord,1) = xyzhex(xcoord,v1l) - xyzavg(xcoord)
311  xyznodes(xcoord,2) = xyzhex(xcoord,v2l) - xyzavg(xcoord)
312  xyznodes(xcoord,3) = xyzhex(xcoord,v3l) - xyzavg(xcoord)
313  xyznodes(xcoord,4) = xyzhex(xcoord,v4l) - xyzavg(xcoord)
314 
315  xyznodes(ycoord,1) = xyzhex(ycoord,v1l) - xyzavg(ycoord)
316  xyznodes(ycoord,2) = xyzhex(ycoord,v2l) - xyzavg(ycoord)
317  xyznodes(ycoord,3) = xyzhex(ycoord,v3l) - xyzavg(ycoord)
318  xyznodes(ycoord,4) = xyzhex(ycoord,v4l) - xyzavg(ycoord)
319 
320  xyznodes(zcoord,1) = xyzhex(zcoord,v1l) - xyzavg(zcoord)
321  xyznodes(zcoord,2) = xyzhex(zcoord,v2l) - xyzavg(zcoord)
322  xyznodes(zcoord,3) = xyzhex(zcoord,v3l) - xyzavg(zcoord)
323  xyznodes(zcoord,4) = xyzhex(zcoord,v4l) - xyzavg(zcoord)
324 
325  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4),fnx,fny,fnz)
326  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4),fcx,fcy,fcz)
327 
328  IF ( global%solverType == solv_implicit_nk ) THEN
329  vsnew = vsnew + (fcx*fnx + fcy*fny + fcz*fnz)
330  ELSE
331  pgrid%gs(ifc) = pgrid%gs(ifc) + (fcx*fnx + fcy*fny + fcz*fnz)
332  END IF ! global%solverType
333 
334  END DO ! ifcl
335  END IF ! v4
336 
337 ! ******************************************************************************
338 ! Interior faces with Old2 grid for implicit code
339 ! ******************************************************************************
340 
341  IF ( global%solverType == solv_implicit_nk ) THEN
342 
343 ! ==============================================================================
344 ! Triangular face, sweeps out prism during grid motion
345 ! ==============================================================================
346 
347  IF ( v4 == vert_none ) THEN ! triangular face
348  xyzpri(xcoord,1) = pxyzold2(xcoord,v1)
349  xyzpri(xcoord,2) = pxyzold2(xcoord,v2)
350  xyzpri(xcoord,3) = pxyzold2(xcoord,v3)
351 
352  xyzpri(ycoord,1) = pxyzold2(ycoord,v1)
353  xyzpri(ycoord,2) = pxyzold2(ycoord,v2)
354  xyzpri(ycoord,3) = pxyzold2(ycoord,v3)
355 
356  xyzpri(zcoord,1) = pxyzold2(zcoord,v1)
357  xyzpri(zcoord,2) = pxyzold2(zcoord,v2)
358  xyzpri(zcoord,3) = pxyzold2(zcoord,v3)
359 
360  xyzpri(xcoord,4) = pxyzold(xcoord,v1)
361  xyzpri(xcoord,5) = pxyzold(xcoord,v2)
362  xyzpri(xcoord,6) = pxyzold(xcoord,v3)
363 
364  xyzpri(ycoord,4) = pxyzold(ycoord,v1)
365  xyzpri(ycoord,5) = pxyzold(ycoord,v2)
366  xyzpri(ycoord,6) = pxyzold(ycoord,v3)
367 
368  xyzpri(zcoord,4) = pxyzold(zcoord,v1)
369  xyzpri(zcoord,5) = pxyzold(zcoord,v2)
370  xyzpri(zcoord,6) = pxyzold(zcoord,v3)
371 
372 ! ------------------------------------------------------------------------------
373 ! Compute average coordinates (approximate centroid)
374 ! ------------------------------------------------------------------------------
375 
376  xyzavg(xcoord) = 0.0_rfreal
377  xyzavg(ycoord) = 0.0_rfreal
378  xyzavg(zcoord) = 0.0_rfreal
379 
380  term = 1.0_rfreal/6.0_rfreal
381 
382  DO iv = 1,6
383  xyzavg(xcoord) = xyzavg(xcoord) + xyzpri(xcoord,iv)
384  xyzavg(ycoord) = xyzavg(ycoord) + xyzpri(ycoord,iv)
385  xyzavg(zcoord) = xyzavg(zcoord) + xyzpri(zcoord,iv)
386  END DO ! iv
387 
388  xyzavg(xcoord) = term*xyzavg(xcoord)
389  xyzavg(ycoord) = term*xyzavg(ycoord)
390  xyzavg(zcoord) = term*xyzavg(zcoord)
391 
392 ! ------------------------------------------------------------------------------
393 ! Loop over the five faces of the prism to compute volume
394 ! ------------------------------------------------------------------------------
395 
396  DO ifcl = 1,5
397  v1l = f2vpri(1,ifcl)
398  v2l = f2vpri(2,ifcl)
399  v3l = f2vpri(3,ifcl)
400  v4l = f2vpri(4,ifcl)
401 
402  IF ( v4l == 0 ) THEN ! triangular face
403  xyznodes(xcoord,1) = xyzpri(xcoord,v1l) - xyzavg(xcoord)
404  xyznodes(xcoord,2) = xyzpri(xcoord,v2l) - xyzavg(xcoord)
405  xyznodes(xcoord,3) = xyzpri(xcoord,v3l) - xyzavg(xcoord)
406 
407  xyznodes(ycoord,1) = xyzpri(ycoord,v1l) - xyzavg(ycoord)
408  xyznodes(ycoord,2) = xyzpri(ycoord,v2l) - xyzavg(ycoord)
409  xyznodes(ycoord,3) = xyzpri(ycoord,v3l) - xyzavg(ycoord)
410 
411  xyznodes(zcoord,1) = xyzpri(zcoord,v1l) - xyzavg(zcoord)
412  xyznodes(zcoord,2) = xyzpri(zcoord,v2l) - xyzavg(zcoord)
413  xyznodes(zcoord,3) = xyzpri(zcoord,v3l) - xyzavg(zcoord)
414 
415  CALL facevectortria(xyznodes(xcoord:zcoord,1:3),fnx,fny,fnz)
416  CALL facecentroidtria(xyznodes(xcoord:zcoord,1:3),fcx,fcy,fcz)
417  ELSE
418  xyznodes(xcoord,1) = xyzpri(xcoord,v1l) - xyzavg(xcoord)
419  xyznodes(xcoord,2) = xyzpri(xcoord,v2l) - xyzavg(xcoord)
420  xyznodes(xcoord,3) = xyzpri(xcoord,v3l) - xyzavg(xcoord)
421  xyznodes(xcoord,4) = xyzpri(xcoord,v4l) - xyzavg(xcoord)
422 
423  xyznodes(ycoord,1) = xyzpri(ycoord,v1l) - xyzavg(ycoord)
424  xyznodes(ycoord,2) = xyzpri(ycoord,v2l) - xyzavg(ycoord)
425  xyznodes(ycoord,3) = xyzpri(ycoord,v3l) - xyzavg(ycoord)
426  xyznodes(ycoord,4) = xyzpri(ycoord,v4l) - xyzavg(ycoord)
427 
428  xyznodes(zcoord,1) = xyzpri(zcoord,v1l) - xyzavg(zcoord)
429  xyznodes(zcoord,2) = xyzpri(zcoord,v2l) - xyzavg(zcoord)
430  xyznodes(zcoord,3) = xyzpri(zcoord,v3l) - xyzavg(zcoord)
431  xyznodes(zcoord,4) = xyzpri(zcoord,v4l) - xyzavg(zcoord)
432 
433  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4),fnx,fny,fnz)
434  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4),fcx,fcy,fcz)
435  END IF ! v4l
436 
437  vsold = vsold + (fcx*fnx + fcy*fny + fcz*fnz)
438 
439  END DO ! ifcl
440 
441 ! ==============================================================================
442 ! Quadrilateral face, sweeps out hexahedron during grid motion
443 ! ==============================================================================
444 
445  ELSE ! quadrilateral face
446  xyzhex(xcoord,1) = pxyzold2(xcoord,v1)
447  xyzhex(xcoord,2) = pxyzold2(xcoord,v2)
448  xyzhex(xcoord,3) = pxyzold2(xcoord,v3)
449  xyzhex(xcoord,4) = pxyzold2(xcoord,v4)
450 
451  xyzhex(ycoord,1) = pxyzold2(ycoord,v1)
452  xyzhex(ycoord,2) = pxyzold2(ycoord,v2)
453  xyzhex(ycoord,3) = pxyzold2(ycoord,v3)
454  xyzhex(ycoord,4) = pxyzold2(ycoord,v4)
455 
456  xyzhex(zcoord,1) = pxyzold2(zcoord,v1)
457  xyzhex(zcoord,2) = pxyzold2(zcoord,v2)
458  xyzhex(zcoord,3) = pxyzold2(zcoord,v3)
459  xyzhex(zcoord,4) = pxyzold2(zcoord,v4)
460 
461  xyzhex(xcoord,5) = pxyzold(xcoord,v1)
462  xyzhex(xcoord,6) = pxyzold(xcoord,v2)
463  xyzhex(xcoord,7) = pxyzold(xcoord,v3)
464  xyzhex(xcoord,8) = pxyzold(xcoord,v4)
465 
466  xyzhex(ycoord,5) = pxyzold(ycoord,v1)
467  xyzhex(ycoord,6) = pxyzold(ycoord,v2)
468  xyzhex(ycoord,7) = pxyzold(ycoord,v3)
469  xyzhex(ycoord,8) = pxyzold(ycoord,v4)
470 
471  xyzhex(zcoord,5) = pxyzold(zcoord,v1)
472  xyzhex(zcoord,6) = pxyzold(zcoord,v2)
473  xyzhex(zcoord,7) = pxyzold(zcoord,v3)
474  xyzhex(zcoord,8) = pxyzold(zcoord,v4)
475 
476 ! ------------------------------------------------------------------------------
477 ! Compute average coordinates (approximate centroid)
478 ! ------------------------------------------------------------------------------
479 
480  xyzavg(xcoord) = 0.0_rfreal
481  xyzavg(ycoord) = 0.0_rfreal
482  xyzavg(zcoord) = 0.0_rfreal
483 
484  term = 1.0_rfreal/8.0_rfreal
485 
486  DO iv = 1,8
487  xyzavg(xcoord) = xyzavg(xcoord) + xyzhex(xcoord,iv)
488  xyzavg(ycoord) = xyzavg(ycoord) + xyzhex(ycoord,iv)
489  xyzavg(zcoord) = xyzavg(zcoord) + xyzhex(zcoord,iv)
490  END DO ! iv
491 
492  xyzavg(xcoord) = term*xyzavg(xcoord)
493  xyzavg(ycoord) = term*xyzavg(ycoord)
494  xyzavg(zcoord) = term*xyzavg(zcoord)
495 
496 ! ------------------------------------------------------------------------------
497 ! Loop over the six faces of the prism to compute volume
498 ! ------------------------------------------------------------------------------
499 
500  DO ifcl = 1,6
501  v1l = f2vhex(1,ifcl)
502  v2l = f2vhex(2,ifcl)
503  v3l = f2vhex(3,ifcl)
504  v4l = f2vhex(4,ifcl)
505 
506  xyznodes(xcoord,1) = xyzhex(xcoord,v1l) - xyzavg(xcoord)
507  xyznodes(xcoord,2) = xyzhex(xcoord,v2l) - xyzavg(xcoord)
508  xyznodes(xcoord,3) = xyzhex(xcoord,v3l) - xyzavg(xcoord)
509  xyznodes(xcoord,4) = xyzhex(xcoord,v4l) - xyzavg(xcoord)
510 
511  xyznodes(ycoord,1) = xyzhex(ycoord,v1l) - xyzavg(ycoord)
512  xyznodes(ycoord,2) = xyzhex(ycoord,v2l) - xyzavg(ycoord)
513  xyznodes(ycoord,3) = xyzhex(ycoord,v3l) - xyzavg(ycoord)
514  xyznodes(ycoord,4) = xyzhex(ycoord,v4l) - xyzavg(ycoord)
515 
516  xyznodes(zcoord,1) = xyzhex(zcoord,v1l) - xyzavg(zcoord)
517  xyznodes(zcoord,2) = xyzhex(zcoord,v2l) - xyzavg(zcoord)
518  xyznodes(zcoord,3) = xyzhex(zcoord,v3l) - xyzavg(zcoord)
519  xyznodes(zcoord,4) = xyzhex(zcoord,v4l) - xyzavg(zcoord)
520 
521  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4),fnx,fny,fnz)
522  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4),fcx,fcy,fcz)
523 
524  vsold = vsold + (fcx*fnx + fcy*fny + fcz*fnz)
525  END DO ! ifcl
526  END IF ! v4
527  END IF ! global%solverType
528 
529 ! ==============================================================================
530 ! Normalize grid speed
531 ! ==============================================================================
532 
533  IF ( global%solverType == solv_implicit_nk ) THEN
534  vsnew = thrd*vsnew
535  vsold = thrd*vsold
536 
537  pgrid%gs(ifc) = (3.0_rfreal*vsnew - vsold)/ &
538  (2.0_rfreal*pgrid%fn(xyzmag,ifc)*global%dtMin)
539  ELSE
540  pgrid%gs(ifc) = thrd*pgrid%gs(ifc)/(pgrid%fn(xyzmag,ifc)*global%dtMin)
541  END IF ! global%solverType
542  END DO ! ifc
543 
544 ! ******************************************************************************
545 ! Write out extrema for interior faces
546 ! ******************************************************************************
547 
548  IF ( global%myProcid == masterproc .AND. &
549  global%verbLevel >= verbose_high ) THEN
550  IF ( pgrid%nFaces > 0 ) THEN ! Have gm test case with 0 faces...
551  WRITE(stdout,'(A,5X,A,2(1X,E24.15),2(1X,I9))') &
552  solver_name,'Interior:', &
553  minval(pgrid%gs),maxval(pgrid%gs), &
554  minloc(pgrid%gs),maxloc(pgrid%gs)
555  END IF ! pGrid%nFacesTot
556  END IF ! global%myProcid
557 
558 ! ******************************************************************************
559 ! Boundary faces
560 ! ******************************************************************************
561 
562  DO ipatch = 1,pgrid%nPatches
563  ppatch => pregion%patches(ipatch)
564 
565 ! ==============================================================================
566 ! Moving faces
567 ! ==============================================================================
568 
569  IF ( ppatch%movePatchDir /= movepatch_dir_none ) THEN
570  DO ifc = 1,ppatch%nBFaces
571  v1 = ppatch%bv(ppatch%bf2v(1,ifc))
572  v2 = ppatch%bv(ppatch%bf2v(2,ifc))
573  v3 = ppatch%bv(ppatch%bf2v(3,ifc))
574 
575  ppatch%gs(ifc) = 0.0_rfreal
576 
577  IF ( global%solverType == solv_implicit_nk ) THEN
578  vsnew = 0.0_rfreal
579  vsold = 0.0_rfreal
580  END IF ! global%solverType
581 
582 ! ------------------------------------------------------------------------------
583 ! Triangular face, sweeps out prism during grid motion
584 ! ------------------------------------------------------------------------------
585 
586  IF ( ppatch%bf2v(4,ifc) == vert_none ) THEN ! triangular face
587  xyzpri(xcoord,1) = pxyzold(xcoord,v1)
588  xyzpri(xcoord,2) = pxyzold(xcoord,v2)
589  xyzpri(xcoord,3) = pxyzold(xcoord,v3)
590 
591  xyzpri(ycoord,1) = pxyzold(ycoord,v1)
592  xyzpri(ycoord,2) = pxyzold(ycoord,v2)
593  xyzpri(ycoord,3) = pxyzold(ycoord,v3)
594 
595  xyzpri(zcoord,1) = pxyzold(zcoord,v1)
596  xyzpri(zcoord,2) = pxyzold(zcoord,v2)
597  xyzpri(zcoord,3) = pxyzold(zcoord,v3)
598 
599  xyzpri(xcoord,4) = pxyz(xcoord,v1)
600  xyzpri(xcoord,5) = pxyz(xcoord,v2)
601  xyzpri(xcoord,6) = pxyz(xcoord,v3)
602 
603  xyzpri(ycoord,4) = pxyz(ycoord,v1)
604  xyzpri(ycoord,5) = pxyz(ycoord,v2)
605  xyzpri(ycoord,6) = pxyz(ycoord,v3)
606 
607  xyzpri(zcoord,4) = pxyz(zcoord,v1)
608  xyzpri(zcoord,5) = pxyz(zcoord,v2)
609  xyzpri(zcoord,6) = pxyz(zcoord,v3)
610 
611 ! ------- Compute average coordinates (approximate centroid)
612 
613  xyzavg(xcoord) = 0.0_rfreal
614  xyzavg(ycoord) = 0.0_rfreal
615  xyzavg(zcoord) = 0.0_rfreal
616 
617  term = 1.0_rfreal/6.0_rfreal
618 
619  DO iv = 1,6
620  xyzavg(xcoord) = xyzavg(xcoord) + xyzpri(xcoord,iv)
621  xyzavg(ycoord) = xyzavg(ycoord) + xyzpri(ycoord,iv)
622  xyzavg(zcoord) = xyzavg(zcoord) + xyzpri(zcoord,iv)
623  END DO ! iv
624 
625  xyzavg(xcoord) = term*xyzavg(xcoord)
626  xyzavg(ycoord) = term*xyzavg(ycoord)
627  xyzavg(zcoord) = term*xyzavg(zcoord)
628 
629 ! ------- Loop over the five faces of the prism to compute volume
630 
631  DO ifcl = 1,5
632  v1l = f2vpri(1,ifcl)
633  v2l = f2vpri(2,ifcl)
634  v3l = f2vpri(3,ifcl)
635  v4l = f2vpri(4,ifcl)
636 
637  IF ( v4l == 0 ) THEN ! triangular face
638  xyznodes(xcoord,1) = xyzpri(xcoord,v1l) - xyzavg(xcoord)
639  xyznodes(xcoord,2) = xyzpri(xcoord,v2l) - xyzavg(xcoord)
640  xyznodes(xcoord,3) = xyzpri(xcoord,v3l) - xyzavg(xcoord)
641 
642  xyznodes(ycoord,1) = xyzpri(ycoord,v1l) - xyzavg(ycoord)
643  xyznodes(ycoord,2) = xyzpri(ycoord,v2l) - xyzavg(ycoord)
644  xyznodes(ycoord,3) = xyzpri(ycoord,v3l) - xyzavg(ycoord)
645 
646  xyznodes(zcoord,1) = xyzpri(zcoord,v1l) - xyzavg(zcoord)
647  xyznodes(zcoord,2) = xyzpri(zcoord,v2l) - xyzavg(zcoord)
648  xyznodes(zcoord,3) = xyzpri(zcoord,v3l) - xyzavg(zcoord)
649 
650  CALL facevectortria(xyznodes(xcoord:zcoord,1:3),fnx,fny,fnz)
651  CALL facecentroidtria(xyznodes(xcoord:zcoord,1:3),fcx,fcy,fcz)
652  ELSE
653  xyznodes(xcoord,1) = xyzpri(xcoord,v1l) - xyzavg(xcoord)
654  xyznodes(xcoord,2) = xyzpri(xcoord,v2l) - xyzavg(xcoord)
655  xyznodes(xcoord,3) = xyzpri(xcoord,v3l) - xyzavg(xcoord)
656  xyznodes(xcoord,4) = xyzpri(xcoord,v4l) - xyzavg(xcoord)
657 
658  xyznodes(ycoord,1) = xyzpri(ycoord,v1l) - xyzavg(ycoord)
659  xyznodes(ycoord,2) = xyzpri(ycoord,v2l) - xyzavg(ycoord)
660  xyznodes(ycoord,3) = xyzpri(ycoord,v3l) - xyzavg(ycoord)
661  xyznodes(ycoord,4) = xyzpri(ycoord,v4l) - xyzavg(ycoord)
662 
663  xyznodes(zcoord,1) = xyzpri(zcoord,v1l) - xyzavg(zcoord)
664  xyznodes(zcoord,2) = xyzpri(zcoord,v2l) - xyzavg(zcoord)
665  xyznodes(zcoord,3) = xyzpri(zcoord,v3l) - xyzavg(zcoord)
666  xyznodes(zcoord,4) = xyzpri(zcoord,v4l) - xyzavg(zcoord)
667 
668  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4),fnx,fny,fnz)
669  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4),fcx,fcy,fcz)
670  END IF ! v4l
671 
672  IF ( global%solverType == solv_implicit_nk ) THEN
673  vsnew = vsnew + (fcx*fnx + fcy*fny + fcz*fnz)
674  ELSE
675  ppatch%gs(ifc) = ppatch%gs(ifc) + (fcx*fnx + fcy*fny + fcz*fnz)
676  END IF ! global%solverType
677 
678  END DO ! ifcl
679 
680 ! ------------------------------------------------------------------------------
681 ! Quadrilateral face, sweeps out hexahedron during grid motion
682 ! ------------------------------------------------------------------------------
683 
684  ELSE ! quadrilateral face
685  v4 = ppatch%bv(ppatch%bf2v(4,ifc))
686 
687  xyzhex(xcoord,1) = pxyzold(xcoord,v1)
688  xyzhex(xcoord,2) = pxyzold(xcoord,v2)
689  xyzhex(xcoord,3) = pxyzold(xcoord,v3)
690  xyzhex(xcoord,4) = pxyzold(xcoord,v4)
691 
692  xyzhex(ycoord,1) = pxyzold(ycoord,v1)
693  xyzhex(ycoord,2) = pxyzold(ycoord,v2)
694  xyzhex(ycoord,3) = pxyzold(ycoord,v3)
695  xyzhex(ycoord,4) = pxyzold(ycoord,v4)
696 
697  xyzhex(zcoord,1) = pxyzold(zcoord,v1)
698  xyzhex(zcoord,2) = pxyzold(zcoord,v2)
699  xyzhex(zcoord,3) = pxyzold(zcoord,v3)
700  xyzhex(zcoord,4) = pxyzold(zcoord,v4)
701 
702  xyzhex(xcoord,5) = pxyz(xcoord,v1)
703  xyzhex(xcoord,6) = pxyz(xcoord,v2)
704  xyzhex(xcoord,7) = pxyz(xcoord,v3)
705  xyzhex(xcoord,8) = pxyz(xcoord,v4)
706 
707  xyzhex(ycoord,5) = pxyz(ycoord,v1)
708  xyzhex(ycoord,6) = pxyz(ycoord,v2)
709  xyzhex(ycoord,7) = pxyz(ycoord,v3)
710  xyzhex(ycoord,8) = pxyz(ycoord,v4)
711 
712  xyzhex(zcoord,5) = pxyz(zcoord,v1)
713  xyzhex(zcoord,6) = pxyz(zcoord,v2)
714  xyzhex(zcoord,7) = pxyz(zcoord,v3)
715  xyzhex(zcoord,8) = pxyz(zcoord,v4)
716 
717 ! ------- Compute average coordinates (approximate centroid)
718 
719  xyzavg(xcoord) = 0.0_rfreal
720  xyzavg(ycoord) = 0.0_rfreal
721  xyzavg(zcoord) = 0.0_rfreal
722 
723  term = 1.0_rfreal/8.0_rfreal
724 
725  DO iv = 1,8
726  xyzavg(xcoord) = xyzavg(xcoord) + xyzhex(xcoord,iv)
727  xyzavg(ycoord) = xyzavg(ycoord) + xyzhex(ycoord,iv)
728  xyzavg(zcoord) = xyzavg(zcoord) + xyzhex(zcoord,iv)
729  END DO ! iv
730 
731  xyzavg(xcoord) = term*xyzavg(xcoord)
732  xyzavg(ycoord) = term*xyzavg(ycoord)
733  xyzavg(zcoord) = term*xyzavg(zcoord)
734 
735 !-------- Loop over the six faces of the hexahedron to compute volume
736 
737  DO ifcl = 1,6
738  v1l = f2vhex(1,ifcl)
739  v2l = f2vhex(2,ifcl)
740  v3l = f2vhex(3,ifcl)
741  v4l = f2vhex(4,ifcl)
742 
743  xyznodes(xcoord,1) = xyzhex(xcoord,v1l) - xyzavg(xcoord)
744  xyznodes(xcoord,2) = xyzhex(xcoord,v2l) - xyzavg(xcoord)
745  xyznodes(xcoord,3) = xyzhex(xcoord,v3l) - xyzavg(xcoord)
746  xyznodes(xcoord,4) = xyzhex(xcoord,v4l) - xyzavg(xcoord)
747 
748  xyznodes(ycoord,1) = xyzhex(ycoord,v1l) - xyzavg(ycoord)
749  xyznodes(ycoord,2) = xyzhex(ycoord,v2l) - xyzavg(ycoord)
750  xyznodes(ycoord,3) = xyzhex(ycoord,v3l) - xyzavg(ycoord)
751  xyznodes(ycoord,4) = xyzhex(ycoord,v4l) - xyzavg(ycoord)
752 
753  xyznodes(zcoord,1) = xyzhex(zcoord,v1l) - xyzavg(zcoord)
754  xyznodes(zcoord,2) = xyzhex(zcoord,v2l) - xyzavg(zcoord)
755  xyznodes(zcoord,3) = xyzhex(zcoord,v3l) - xyzavg(zcoord)
756  xyznodes(zcoord,4) = xyzhex(zcoord,v4l) - xyzavg(zcoord)
757 
758  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4),fnx,fny,fnz)
759  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4),fcx,fcy,fcz)
760 
761  IF ( global%solverType == solv_implicit_nk ) THEN
762  vsnew = vsnew + (fcx*fnx + fcy*fny + fcz*fnz)
763  ELSE
764  ppatch%gs(ifc) = ppatch%gs(ifc) + (fcx*fnx + fcy*fny + fcz*fnz)
765  END IF ! global%solverType
766 
767  END DO ! ifcl
768  END IF ! v4
769 
770 ! ******************************************************************************
771 ! Boundary faces with Old2 grid for implicit code
772 ! ******************************************************************************
773 
774  IF ( global%solverType == solv_implicit_nk ) THEN
775 
776 ! ------------------------------------------------------------------------------
777 ! Triangular face, sweeps out prism during grid motion
778 ! ------------------------------------------------------------------------------
779 
780  IF ( ppatch%bf2v(4,ifc) == vert_none ) THEN ! triangular face
781  xyzpri(xcoord,1) = pxyzold2(xcoord,v1)
782  xyzpri(xcoord,2) = pxyzold2(xcoord,v2)
783  xyzpri(xcoord,3) = pxyzold2(xcoord,v3)
784 
785  xyzpri(ycoord,1) = pxyzold2(ycoord,v1)
786  xyzpri(ycoord,2) = pxyzold2(ycoord,v2)
787  xyzpri(ycoord,3) = pxyzold2(ycoord,v3)
788 
789  xyzpri(zcoord,1) = pxyzold2(zcoord,v1)
790  xyzpri(zcoord,2) = pxyzold2(zcoord,v2)
791  xyzpri(zcoord,3) = pxyzold2(zcoord,v3)
792 
793  xyzpri(xcoord,4) = pxyzold(xcoord,v1)
794  xyzpri(xcoord,5) = pxyzold(xcoord,v2)
795  xyzpri(xcoord,6) = pxyzold(xcoord,v3)
796 
797  xyzpri(ycoord,4) = pxyzold(ycoord,v1)
798  xyzpri(ycoord,5) = pxyzold(ycoord,v2)
799  xyzpri(ycoord,6) = pxyzold(ycoord,v3)
800 
801  xyzpri(zcoord,4) = pxyzold(zcoord,v1)
802  xyzpri(zcoord,5) = pxyzold(zcoord,v2)
803  xyzpri(zcoord,6) = pxyzold(zcoord,v3)
804 
805 ! ------- Compute average coordinates (approximate centroid)
806 
807  xyzavg(xcoord) = 0.0_rfreal
808  xyzavg(ycoord) = 0.0_rfreal
809  xyzavg(zcoord) = 0.0_rfreal
810 
811  term = 1.0_rfreal/6.0_rfreal
812 
813  DO iv = 1,6
814  xyzavg(xcoord) = xyzavg(xcoord) + xyzpri(xcoord,iv)
815  xyzavg(ycoord) = xyzavg(ycoord) + xyzpri(ycoord,iv)
816  xyzavg(zcoord) = xyzavg(zcoord) + xyzpri(zcoord,iv)
817  END DO ! iv
818 
819  xyzavg(xcoord) = term*xyzavg(xcoord)
820  xyzavg(ycoord) = term*xyzavg(ycoord)
821  xyzavg(zcoord) = term*xyzavg(zcoord)
822 
823 ! ------- Loop over the five faces of the prism to compute volume
824 
825  DO ifcl = 1,5
826  v1l = f2vpri(1,ifcl)
827  v2l = f2vpri(2,ifcl)
828  v3l = f2vpri(3,ifcl)
829  v4l = f2vpri(4,ifcl)
830 
831  IF ( v4l == 0 ) THEN ! triangular face
832  xyznodes(xcoord,1) = xyzpri(xcoord,v1l) - xyzavg(xcoord)
833  xyznodes(xcoord,2) = xyzpri(xcoord,v2l) - xyzavg(xcoord)
834  xyznodes(xcoord,3) = xyzpri(xcoord,v3l) - xyzavg(xcoord)
835 
836  xyznodes(ycoord,1) = xyzpri(ycoord,v1l) - xyzavg(ycoord)
837  xyznodes(ycoord,2) = xyzpri(ycoord,v2l) - xyzavg(ycoord)
838  xyznodes(ycoord,3) = xyzpri(ycoord,v3l) - xyzavg(ycoord)
839 
840  xyznodes(zcoord,1) = xyzpri(zcoord,v1l) - xyzavg(zcoord)
841  xyznodes(zcoord,2) = xyzpri(zcoord,v2l) - xyzavg(zcoord)
842  xyznodes(zcoord,3) = xyzpri(zcoord,v3l) - xyzavg(zcoord)
843 
844  CALL facevectortria(xyznodes(xcoord:zcoord,1:3),fnx,fny,fnz)
845  CALL facecentroidtria(xyznodes(xcoord:zcoord,1:3),fcx,fcy,fcz)
846  ELSE
847  xyznodes(xcoord,1) = xyzpri(xcoord,v1l) - xyzavg(xcoord)
848  xyznodes(xcoord,2) = xyzpri(xcoord,v2l) - xyzavg(xcoord)
849  xyznodes(xcoord,3) = xyzpri(xcoord,v3l) - xyzavg(xcoord)
850  xyznodes(xcoord,4) = xyzpri(xcoord,v4l) - xyzavg(xcoord)
851 
852  xyznodes(ycoord,1) = xyzpri(ycoord,v1l) - xyzavg(ycoord)
853  xyznodes(ycoord,2) = xyzpri(ycoord,v2l) - xyzavg(ycoord)
854  xyznodes(ycoord,3) = xyzpri(ycoord,v3l) - xyzavg(ycoord)
855  xyznodes(ycoord,4) = xyzpri(ycoord,v4l) - xyzavg(ycoord)
856 
857  xyznodes(zcoord,1) = xyzpri(zcoord,v1l) - xyzavg(zcoord)
858  xyznodes(zcoord,2) = xyzpri(zcoord,v2l) - xyzavg(zcoord)
859  xyznodes(zcoord,3) = xyzpri(zcoord,v3l) - xyzavg(zcoord)
860  xyznodes(zcoord,4) = xyzpri(zcoord,v4l) - xyzavg(zcoord)
861 
862  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4),fnx,fny,fnz)
863  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4),fcx,fcy,fcz)
864  END IF ! v4l
865 
866  vsold = vsold + (fcx*fnx + fcy*fny + fcz*fnz)
867 
868  END DO ! ifcl
869 
870 ! ------------------------------------------------------------------------------
871 ! Quadrilateral face, sweeps out hexahedron during grid motion
872 ! ------------------------------------------------------------------------------
873 
874  ELSE ! quadrilateral face
875  v4 = ppatch%bv(ppatch%bf2v(4,ifc))
876 
877  xyzhex(xcoord,1) = pxyzold(xcoord,v1)
878  xyzhex(xcoord,2) = pxyzold(xcoord,v2)
879  xyzhex(xcoord,3) = pxyzold(xcoord,v3)
880  xyzhex(xcoord,4) = pxyzold(xcoord,v4)
881 
882  xyzhex(ycoord,1) = pxyzold(ycoord,v1)
883  xyzhex(ycoord,2) = pxyzold(ycoord,v2)
884  xyzhex(ycoord,3) = pxyzold(ycoord,v3)
885  xyzhex(ycoord,4) = pxyzold(ycoord,v4)
886 
887  xyzhex(zcoord,1) = pxyzold(zcoord,v1)
888  xyzhex(zcoord,2) = pxyzold(zcoord,v2)
889  xyzhex(zcoord,3) = pxyzold(zcoord,v3)
890  xyzhex(zcoord,4) = pxyzold(zcoord,v4)
891 
892  xyzhex(xcoord,5) = pxyz(xcoord,v1)
893  xyzhex(xcoord,6) = pxyz(xcoord,v2)
894  xyzhex(xcoord,7) = pxyz(xcoord,v3)
895  xyzhex(xcoord,8) = pxyz(xcoord,v4)
896 
897  xyzhex(ycoord,5) = pxyz(ycoord,v1)
898  xyzhex(ycoord,6) = pxyz(ycoord,v2)
899  xyzhex(ycoord,7) = pxyz(ycoord,v3)
900  xyzhex(ycoord,8) = pxyz(ycoord,v4)
901 
902  xyzhex(zcoord,5) = pxyz(zcoord,v1)
903  xyzhex(zcoord,6) = pxyz(zcoord,v2)
904  xyzhex(zcoord,7) = pxyz(zcoord,v3)
905  xyzhex(zcoord,8) = pxyz(zcoord,v4)
906 
907 ! ------- Compute average coordinates (approximate centroid)
908 
909  xyzavg(xcoord) = 0.0_rfreal
910  xyzavg(ycoord) = 0.0_rfreal
911  xyzavg(zcoord) = 0.0_rfreal
912 
913  term = 1.0_rfreal/8.0_rfreal
914 
915  DO iv = 1,8
916  xyzavg(xcoord) = xyzavg(xcoord) + xyzhex(xcoord,iv)
917  xyzavg(ycoord) = xyzavg(ycoord) + xyzhex(ycoord,iv)
918  xyzavg(zcoord) = xyzavg(zcoord) + xyzhex(zcoord,iv)
919  END DO ! iv
920 
921  xyzavg(xcoord) = term*xyzavg(xcoord)
922  xyzavg(ycoord) = term*xyzavg(ycoord)
923  xyzavg(zcoord) = term*xyzavg(zcoord)
924 
925 !-------- Loop over the six faces of the hexahedron to compute volume
926 
927  DO ifcl = 1,6
928  v1l = f2vhex(1,ifcl)
929  v2l = f2vhex(2,ifcl)
930  v3l = f2vhex(3,ifcl)
931  v4l = f2vhex(4,ifcl)
932 
933  xyznodes(xcoord,1) = xyzhex(xcoord,v1l) - xyzavg(xcoord)
934  xyznodes(xcoord,2) = xyzhex(xcoord,v2l) - xyzavg(xcoord)
935  xyznodes(xcoord,3) = xyzhex(xcoord,v3l) - xyzavg(xcoord)
936  xyznodes(xcoord,4) = xyzhex(xcoord,v4l) - xyzavg(xcoord)
937 
938  xyznodes(ycoord,1) = xyzhex(ycoord,v1l) - xyzavg(ycoord)
939  xyznodes(ycoord,2) = xyzhex(ycoord,v2l) - xyzavg(ycoord)
940  xyznodes(ycoord,3) = xyzhex(ycoord,v3l) - xyzavg(ycoord)
941  xyznodes(ycoord,4) = xyzhex(ycoord,v4l) - xyzavg(ycoord)
942 
943  xyznodes(zcoord,1) = xyzhex(zcoord,v1l) - xyzavg(zcoord)
944  xyznodes(zcoord,2) = xyzhex(zcoord,v2l) - xyzavg(zcoord)
945  xyznodes(zcoord,3) = xyzhex(zcoord,v3l) - xyzavg(zcoord)
946  xyznodes(zcoord,4) = xyzhex(zcoord,v4l) - xyzavg(zcoord)
947 
948  CALL facevectorquad(xyznodes(xcoord:zcoord,1:4),fnx,fny,fnz)
949  CALL facecentroidquad(xyznodes(xcoord:zcoord,1:4),fcx,fcy,fcz)
950 
951  vsold = vsold + (fcx*fnx + fcy*fny + fcz*fnz)
952 
953  END DO ! ifcl
954  END IF ! v4
955 
956  END IF ! global%solverType
957 
958 ! ------------------------------------------------------------------------------
959 ! Normalize grid speed
960 ! ------------------------------------------------------------------------------
961 
962  IF ( global%solverType == solv_implicit_nk ) THEN
963  vsnew = thrd*vsnew
964  vsold = thrd*vsold
965  ppatch%gs(ifc) = (3.0_rfreal*vsnew - vsold)/(2.0_rfreal*ppatch%fn(xyzmag,ifc)*global%dtMin)
966  ELSE
967  ppatch%gs(ifc) = thrd*ppatch%gs(ifc)/ (ppatch%fn(xyzmag,ifc)*global%dtMin)
968  END IF ! global%solverType
969 
970  END DO ! ifc
971 
972 ! ==============================================================================
973 ! Non-moving patches (should have precisely zero grid speed, so enforce)
974 ! ==============================================================================
975 
976  ELSE
977  DO ifc = 1,ppatch%nBFaces
978  ppatch%gs(ifc) = 0.0_rfreal
979  END DO ! ifc
980  END IF ! pPatch%movePatchDir
981 
982 ! ==============================================================================
983 ! Write out extrema for patches with actual faces
984 ! ==============================================================================
985 
986  IF ( global%myProcid == masterproc .AND. &
987  global%verbLevel >= verbose_high ) THEN
988  IF ( ppatch%nBFaces > 0 ) THEN ! can have zero faces
989  WRITE(stdout,'(A,5X,A,I3,A,2(1X,E24.15),2(1X,I9))') &
990  solver_name,'Patch',ipatch,':', &
991  minval(ppatch%gs(1:ppatch%nBFaces)), &
992  maxval(ppatch%gs(1:ppatch%nBFaces)), &
993  minloc(ppatch%gs(1:ppatch%nBFaces)), &
994  maxloc(ppatch%gs(1:ppatch%nBFaces))
995  END IF ! pPatch%nBFaces
996  END IF ! global%myProcid
997  END DO ! iPatch
998 
999 ! ******************************************************************************
1000 ! End
1001 ! ******************************************************************************
1002 
1003  IF ( global%myProcid == masterproc .AND. &
1004  global%verbLevel >= verbose_high ) THEN
1005  WRITE(stdout,'(A,1X,A)') solver_name,'Computing grid speeds done.'
1006  END IF ! global%myProcid
1007 
1008  CALL deregisterfunction(global)
1009 
1010 
1011 END SUBROUTINE rflu_computegridspeeds
1012 
1013 ! ******************************************************************************
1014 !
1015 ! RCS Revision history:
1016 !
1017 ! $Log: RFLU_ComputeGridSpeeds.F90,v $
1018 ! Revision 1.9 2008/12/06 08:44:29 mtcampbe
1019 ! Updated license.
1020 !
1021 ! Revision 1.8 2008/11/19 22:17:42 mtcampbe
1022 ! Added Illinois Open Source License/Copyright
1023 !
1024 ! Revision 1.7 2006/02/08 21:31:48 hdewey2
1025 ! Added grid speed computation for implicit solver
1026 !
1027 ! Revision 1.6 2005/06/09 20:22:58 haselbac
1028 ! Replaced movePatch by movePatchDir
1029 !
1030 ! Revision 1.5 2004/10/19 19:41:18 haselbac
1031 ! Cosmetics only
1032 !
1033 ! Revision 1.4 2003/03/15 18:28:01 haselbac
1034 ! Changed loop limits, added VERT_NONE, other changes for || gm
1035 !
1036 ! Revision 1.3 2003/01/28 14:31:27 haselbac
1037 ! Set gs to zero on non-moving boundaries, use scaled coordinates, clean-up
1038 !
1039 ! Revision 1.2 2002/11/08 21:29:45 haselbac
1040 ! Some cosmetics and clean-up
1041 !
1042 ! Revision 1.1 2002/10/27 19:20:08 haselbac
1043 ! Initial revision
1044 !
1045 ! ******************************************************************************
1046 
1047 
1048 
1049 
1050 
1051 
1052 
rflu_modgrid
Definition: RFLU_ModGrid.F90:39
facevectortria
subroutine facevectortria(xyzNodes, fVecX, fVecY, fVecZ)
Definition: FaceVector.F90:49
Kdtree_d::ExtPoint::term
void term()
Definition: kdtree_d.h:237
moderror::registerfunction
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
modgrid::t_grid
Definition: ModGrid.F90:56
rflu_computegridspeeds
subroutine rflu_computegridspeeds(pRegion)
Definition: RFLU_ComputeGridSpeeds.F90:55
modglobal::t_global
Definition: ModGlobal.F90:50
facevectorquad
subroutine facevectorquad(xyzNodes, fVecX, fVecY, fVecZ)
Definition: FaceVector.F90:80
moderror::deregisterfunction
subroutine deregisterfunction(global)
Definition: ModError.F90:469
facecentroidtria
subroutine facecentroidtria(xyz, fCenX, fCenY, fCenZ)
Definition: FaceCentroid.F90:47
facecentroidquad
subroutine facecentroidquad(xyz, fCenX, fCenY, fCenZ)
Definition: FaceCentroid.F90:73
modbndpatch::t_patch
Definition: ModBndPatch.F90:91