Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLU_ModDifferentiationCells.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 to differentiate function values at cell centroids.
26 !
27 ! Description: None.
28 !
29 ! Notes: None.
30 !
31 ! ******************************************************************************
32 !
33 ! $Id: RFLU_ModDifferentiationCells.F90,v 1.12 2008/12/06 08:44:21 mtcampbe Exp $
34 !
35 ! Copyright: (c) 2005-2006 by the University of Illinois
36 !
37 ! ******************************************************************************
38 
40 
41  USE modglobal, ONLY: t_global
42  USE modparameters
43  USE moddatatypes
44  USE moderror
45  USE moddatastruct, ONLY: t_region
46  USE modgrid, ONLY: t_grid
47  USE modbndpatch, ONLY: t_patch
48  USE modmpi
49 
51 
52  IMPLICIT NONE
53 
54  PRIVATE
56 
57 ! ******************************************************************************
58 ! Declarations and definitions
59 ! ******************************************************************************
60 
61  CHARACTER(CHRLEN) :: RCSIdentString = &
62  '$RCSfile: RFLU_ModDifferentiationCells.F90,v $ $Revision: 1.12 $'
63 
64 ! ******************************************************************************
65 ! Routines
66 ! ******************************************************************************
67 
68  CONTAINS
69 
70 
71 
72 
73 
74 ! ******************************************************************************
75 !
76 ! Purpose: Compute 1D gradients of any vector or scalar at cell centers.
77 !
78 ! Description: None.
79 !
80 ! Input:
81 ! pRegion Pointer to region data
82 ! iBegVar Beginning index of data in var
83 ! iEndVar Ending index of data in var
84 ! iBegGrad Beginning index of data in grad
85 ! iEndGrad Ending index of data in grad
86 ! var Variables of which gradients are to be determined
87 !
88 ! Output:
89 ! grad Gradients of variables at cell centers
90 !
91 ! Notes: None.
92 !
93 ! ******************************************************************************
94 
95  SUBROUTINE rflu_computegradcells_1d(pRegion,iBegVar,iEndVar,iBegGrad, &
96  iendgrad,var,grad)
97 
99 
100  IMPLICIT NONE
101 
102 ! ******************************************************************************
103 ! Definitions and declarations
104 ! ******************************************************************************
105 
106 ! ==============================================================================
107 ! Arguments
108 ! ==============================================================================
109 
110  INTEGER, INTENT(IN) :: ibegvar,iendvar,ibeggrad,iendgrad
111  REAL(RFREAL), DIMENSION(:,:), POINTER :: var
112  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
113  TYPE(t_region), POINTER :: pregion
114 
115 ! ==============================================================================
116 ! Locals
117 ! ==============================================================================
118 
119  LOGICAL :: icgincludeflag
120  INTEGER :: errorflag,fndir,fndirend,icg,icg2,igrad,isl,ivar,nmembsmax, &
121  nmembs,order
122  REAL(RFREAL), DIMENSION(:), ALLOCATABLE :: locs,wts
123  TYPE(t_global), POINTER :: global
124  TYPE(t_grid), POINTER :: pgrid
125 
126 ! *****************************************************************************
127 ! Start
128 ! *****************************************************************************
129 
130  global => pregion%global
131 
132  CALL registerfunction(global,'RFLU_ComputeGradCells_1D',&
133  'RFLU_ModDifferentiationCells.F90' )
134 
135 #ifdef ROCPROF
136  CALL fprofiler_begins("RFLU::ComputeGradCells_1D")
137 #endif
138 
139  IF ( (iendvar - ibegvar) /= (iendgrad - ibeggrad) ) THEN
140  CALL errorstop(global,err_grad_mismatch,__line__)
141  END IF ! iEndVar
142 
143 ! ******************************************************************************
144 ! Set pointers and variables
145 ! ******************************************************************************
146 
147  pgrid => pregion%grid
148 
149  nmembsmax = pgrid%c2csInfo%nCellMembsMax
150  order = 1 ! Order of derivative
151 
152  icgincludeflag = .true.
153 
154  SELECT CASE ( pregion%mixtInput%dimens )
155  CASE ( 1 )
156  fndirend = xcoord
157  CASE ( 2 )
158  fndirend = ycoord
159  CASE ( 3 )
160  fndirend = zcoord
161  CASE default
162  CALL errorstop(global,err_reached_default,__line__)
163  END SELECT ! pRegion%mixtInput%dimens
164 
165 ! ******************************************************************************
166 ! Allocate temporary memory
167 ! ******************************************************************************
168 
169  ALLOCATE(wts(0:nmembsmax),stat=errorflag)
170  global%error = errorflag
171  IF ( global%error /= err_none ) THEN
172  CALL errorstop(global,err_allocate,__line__,'wts')
173  END IF ! global%error
174 
175  ALLOCATE(locs(0:nmembsmax),stat=errorflag)
176  global%error = errorflag
177  IF ( global%error /= err_none ) THEN
178  CALL errorstop(global,err_allocate,__line__,'locs')
179  END IF ! global%error
180 
181 ! ******************************************************************************
182 ! Compute gradients
183 ! ******************************************************************************
184 
185 ! ==============================================================================
186 ! Include cell icg in stencil
187 ! ==============================================================================
188 
189  IF ( icgincludeflag .EQV. .true. ) THEN
190  DO icg = 1,pgrid%nCellsTot
191  DO igrad = ibeggrad,iendgrad
192  grad(xcoord,igrad,icg) = 0.0_rfreal
193  grad(ycoord,igrad,icg) = 0.0_rfreal
194  grad(zcoord,igrad,icg) = 0.0_rfreal
195  END DO ! iGrad
196 
197  DO fndir = xcoord,fndirend
198  nmembs = pgrid%c2cs1D(fndir,icg)%nCellMembs
199 
200  locs(0) = pgrid%cofg(fndir,icg)
201 
202  DO isl = 1,nmembs
203  icg2 = pgrid%c2cs1D(fndir,icg)%cellMembs(isl)
204 
205  locs(isl) = pgrid%cofg(fndir,icg2)
206  END DO ! isl
207 
208  CALL rflu_computewtsx2c_1d(global,order,nmembs+1,locs(0:nmembs), &
209  pgrid%cofg(fndir,icg),wts(0:nmembs))
210 
211  igrad = ibeggrad
212 
213  DO ivar = ibegvar,iendvar
214  grad(fndir,igrad,icg) = wts(0)*var(ivar,icg)
215 
216  DO isl = 1,nmembs
217  icg2 = pgrid%c2cs1D(fndir,icg)%cellMembs(isl)
218 
219  grad(fndir,igrad,icg) = grad(fndir,igrad,icg) &
220  + wts(isl)*var(ivar,icg2)
221  END DO ! isl
222 
223  igrad = igrad + 1
224  END DO ! iVar
225  END DO ! fnDir
226  END DO ! icg
227 
228 ! ==============================================================================
229 ! Do not include cell icg in stencil
230 ! ==============================================================================
231 
232  ELSE
233  DO icg = 1,pgrid%nCellsTot
234  DO igrad = ibeggrad,iendgrad
235  grad(xcoord,igrad,icg) = 0.0_rfreal
236  grad(ycoord,igrad,icg) = 0.0_rfreal
237  grad(zcoord,igrad,icg) = 0.0_rfreal
238  END DO ! iGrad
239 
240  DO fndir = xcoord,fndirend
241  nmembs = pgrid%c2cs1D(fndir,icg)%nCellMembs
242 
243  DO isl = 1,nmembs
244  icg2 = pgrid%c2cs1D(fndir,icg)%cellMembs(isl)
245 
246  locs(isl) = pgrid%cofg(fndir,icg2)
247  END DO ! isl
248 
249  CALL rflu_computewtsx2c_1d(global,order,nmembs,locs(1:nmembs), &
250  pgrid%cofg(fndir,icg),wts(1:nmembs))
251 
252  igrad = ibeggrad
253 
254  DO ivar = ibegvar,iendvar
255  DO isl = 1,nmembs
256  icg2 = pgrid%c2cs1D(fndir,icg)%cellMembs(isl)
257 
258  grad(fndir,igrad,icg) = grad(fndir,igrad,icg) &
259  + wts(isl)*var(ivar,icg2)
260  END DO ! isl
261 
262  igrad = igrad + 1
263  END DO ! iVar
264  END DO ! fnDir
265  END DO ! icg
266  END IF ! icgIncludeFlag
267 
268 ! ******************************************************************************
269 ! Deallocate temporary memory
270 ! ******************************************************************************
271 
272  DEALLOCATE(wts,stat=errorflag)
273  global%error = errorflag
274  IF ( global%error /= err_none ) THEN
275  CALL errorstop(global,err_deallocate,__line__,'wts')
276  END IF ! global%error
277 
278  DEALLOCATE(locs,stat=errorflag)
279  global%error = errorflag
280  IF ( global%error /= err_none ) THEN
281  CALL errorstop(global,err_deallocate,__line__,'locs')
282  END IF ! global%error
283 
284 ! ******************************************************************************
285 ! End
286 ! ******************************************************************************
287 
288 #ifdef ROCPROF
289  CALL fprofiler_ends("RFLU::ComputeGradCells_1D")
290 #endif
291 
292  CALL deregisterfunction(global)
293 
294  END SUBROUTINE rflu_computegradcells_1d
295 
296 
297 
298 
299 
300 
301 
302 ! ******************************************************************************
303 !
304 ! Purpose: Compute 2D gradients of any vector or scalar at cell centers.
305 !
306 ! Description: None.
307 !
308 ! Input:
309 ! pRegion Pointer to region data
310 ! iBegVar Beginning index of data in var
311 ! iEndVar Ending index of data in var
312 ! iBegGrad Beginning index of data in grad
313 ! iEndGrad Ending index of data in grad
314 ! var Variables of which gradients are to be determined
315 !
316 ! Output:
317 ! grad Gradients of variables at cell centers
318 !
319 ! Notes:
320 ! 1. Use original form of reconstruction, i.e., weights multiplying variable
321 ! differences, because central location is cell, so have variable located
322 ! there. This is in contrast to face gradients and vertex interpolation,
323 ! where need to use modified form of reconstruction.
324 ! 2. If the weighting is changed from inverse-distance to none, then the
325 ! routine RFLU_ComputeStencilMomentsX in RFLU_ModWeights must also be
326 ! adapted.
327 ! 3. Restricted to linear reconstruction for now.
328 !
329 ! ******************************************************************************
330 
331  SUBROUTINE rflu_computegradcells_2d(pRegion,iBegVar,iEndVar,iBegGrad, &
332  iendgrad,var,grad)
333 
334  IMPLICIT NONE
335 
336 ! ******************************************************************************
337 ! Definitions and declarations
338 ! ******************************************************************************
339 
340 ! ==============================================================================
341 ! Arguments
342 ! ==============================================================================
343 
344  INTEGER, INTENT(IN) :: ibegvar,iendvar,ibeggrad,iendgrad
345  REAL(RFREAL), DIMENSION(:,:), POINTER :: var
346  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
347  TYPE(t_region), POINTER :: pregion
348 
349 ! ==============================================================================
350 ! Locals
351 ! ==============================================================================
352 
353  INTEGER :: errorflag,icg,icg2,igrad,isl,ivar
354  REAL(RFREAL) :: c11,c12,c22,dvar,dx,dy,r11,r12,r22,term,term1,term2,wx,wy
355  REAL(RFREAL) :: cofg(xcoord:ycoord)
356  TYPE(t_global), POINTER :: global
357  TYPE(t_grid), POINTER :: pgrid
358 
359 ! *****************************************************************************
360 ! Start
361 ! *****************************************************************************
362 
363  global => pregion%global
364 
365  CALL registerfunction(global,'RFLU_ComputeGradCells_2D',&
366  'RFLU_ModDifferentiationCells.F90' )
367 
368 #ifdef ROCPROF
369  CALL fprofiler_begins("RFLU::ComputeGradCells_2D")
370 #endif
371 
372  IF ( (iendvar - ibegvar) /= (iendgrad - ibeggrad) ) THEN
373  CALL errorstop(global,err_grad_mismatch,__line__)
374  END IF ! iEndVar
375 
376 ! ******************************************************************************
377 ! Set pointers and variables
378 ! ******************************************************************************
379 
380  pgrid => pregion%grid
381 
382 ! ******************************************************************************
383 ! Loop over cells and compute gradients
384 ! ******************************************************************************
385 
386 ! DEBUG
387 ! DO icg = 1,pGrid%nCellsTot
388 ! DO iVar = iBegVar,iEndVar
389 ! var(iVar,icg) = REAL(4*(iVar-1) ,RFREAL) &
390 ! + REAL(4*(iVar-1)+1,RFREAL)*pGrid%cofg(XCOORD,icg) &
391 ! + REAL(4*(iVar-1)+2,RFREAL)*pGrid%cofg(YCOORD,icg)
392 ! END DO ! iVar
393 ! END DO ! icg
394 ! END DEBUG
395 
396  DO icg = 1,pgrid%nCellsTot
397  r11 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_11)
398  r12 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_12)
399  r22 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_22)
400 
401  c11 = 1.0_rfreal/r11
402  c22 = 1.0_rfreal/r22
403 
404  c12 = -c11*r12
405 
406  cofg(xcoord) = pgrid%cofg(xcoord,icg)
407  cofg(ycoord) = pgrid%cofg(ycoord,icg)
408 
409  DO igrad = ibeggrad,iendgrad
410  grad(xcoord,igrad,icg) = 0.0_rfreal
411  grad(ycoord,igrad,icg) = 0.0_rfreal
412  grad(zcoord,igrad,icg) = 0.0_rfreal ! Set z-component to zero
413  END DO ! iVar
414 
415  DO isl = 1,pgrid%c2cs(icg)%nCellMembs
416  icg2 = pgrid%c2cs(icg)%cellMembs(isl)
417 
418  dx = pgrid%cofg(xcoord,icg2) - cofg(xcoord)
419  dy = pgrid%cofg(ycoord,icg2) - cofg(ycoord)
420 
421  term = 1.0_rfreal/sqrt(dx*dx + dy*dy)
422 
423  dx = term*dx
424  dy = term*dy
425 
426  term1 = c11*c11*( dx)
427  term2 = c22*c22*(dy + c12*dx)
428 
429  wx = term*(term1 + c12*term2)
430  wy = term*( term2)
431 
432  igrad = ibeggrad
433 
434  DO ivar = ibegvar,iendvar
435  dvar = var(ivar,icg2) - var(ivar,icg)
436 
437  grad(xcoord,igrad,icg) = grad(xcoord,igrad,icg) + wx*dvar
438  grad(ycoord,igrad,icg) = grad(ycoord,igrad,icg) + wy*dvar
439 
440  igrad = igrad + 1
441  END DO ! iVar
442  END DO ! isl
443  END DO ! icg
444 
445 ! DEBUG
446 ! DO iGrad = iBegGrad,iEndGrad
447 ! WRITE(*,*) iGrad,MINVAL(grad(XCOORD,iGrad,1:pGrid%nCellsTot)), &
448 ! MINVAL(grad(YCOORD,iGrad,1:pGrid%nCellsTot)), &
449 ! MAXVAL(grad(XCOORD,iGrad,1:pGrid%nCellsTot)), &
450 ! MAXVAL(grad(YCOORD,iGrad,1:pGrid%nCellsTot))
451 ! END DO ! iGrad
452 !
453 ! STOP
454 ! END DEBUG
455 
456 ! ******************************************************************************
457 ! End
458 ! ******************************************************************************
459 
460 #ifdef ROCPROF
461  CALL fprofiler_ends("RFLU::ComputeGradCells_2D")
462 #endif
463 
464  CALL deregisterfunction(global)
465 
466  END SUBROUTINE rflu_computegradcells_2d
467 
468 
469 
470 
471 
472 
473 
474 ! ******************************************************************************
475 !
476 ! Purpose: Compute 3D gradients of any vector or scalar at cell centers.
477 !
478 ! Description: None.
479 !
480 ! Input:
481 ! pRegion Pointer to region data
482 ! iBegVar Beginning index of data in var
483 ! iEndVar Ending index of data in var
484 ! iBegGrad Beginning index of data in grad
485 ! iEndGrad Ending index of data in grad
486 ! var Variables of which gradients are to be determined
487 !
488 ! Output:
489 ! grad Gradients of variables at cell centers
490 !
491 ! Notes:
492 ! 1. Use original form of reconstruction, i.e., weights multiplying variable
493 ! differences, because central location is cell, so have variable located
494 ! there. This is in contrast to face gradients and vertex interpolation,
495 ! where need to use modified form of reconstruction.
496 ! 2. If the weighting is changed from inverse-distance to none, then the
497 ! routine RFLU_ComputeStencilMomentsX in RFLU_ModWeights must also be
498 ! adapted.
499 ! 3. Restricted to linear reconstruction for now.
500 !
501 ! ******************************************************************************
502 
503  SUBROUTINE rflu_computegradcells_3d(pRegion,iBegVar,iEndVar,iBegGrad, &
504  iendgrad,var,grad)
505 
506  IMPLICIT NONE
507 
508 ! ******************************************************************************
509 ! Definitions and declarations
510 ! ******************************************************************************
511 
512 ! ==============================================================================
513 ! Arguments
514 ! ==============================================================================
515 
516  INTEGER, INTENT(IN) :: ibegvar,iendvar,ibeggrad,iendgrad
517  REAL(RFREAL), DIMENSION(:,:), POINTER :: var
518  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
519  TYPE(t_region), POINTER :: pregion
520 
521 ! ==============================================================================
522 ! Locals
523 ! ==============================================================================
524 
525  INTEGER :: errorflag,icg,icg2,igrad,isl,ivar
526  REAL(RFREAL) :: c11,c12,c13,c22,c23,c33,dvar,dx,dy,dz,r11,r12,r13,r22, &
527  r23,r33,term,term1,term2,term3,wx,wy,wz
528  REAL(RFREAL) :: cofg(xcoord:zcoord)
529  TYPE(t_global), POINTER :: global
530  TYPE(t_grid), POINTER :: pgrid
531 
532 ! *****************************************************************************
533 ! Start
534 ! *****************************************************************************
535 
536  global => pregion%global
537 
538  CALL registerfunction(global,'RFLU_ComputeGradCells_3D',&
539  'RFLU_ModDifferentiationCells.F90' )
540 
541 #ifdef ROCPROF
542  CALL fprofiler_begins("RFLU::ComputeGradCells_3D")
543 #endif
544 
545  IF ( (iendvar - ibegvar) /= (iendgrad - ibeggrad) ) THEN
546  CALL errorstop(global,err_grad_mismatch,__line__)
547  END IF ! iEndVar
548 
549 ! ******************************************************************************
550 ! Set pointers and variables
551 ! ******************************************************************************
552 
553  pgrid => pregion%grid
554 
555 ! ******************************************************************************
556 ! Loop over cells and compute gradients
557 ! ******************************************************************************
558 
559 ! DEBUG
560 ! DO icg = 1,pGrid%nCellsTot
561 ! DO iVar = iBegVar,iEndVar
562 ! var(iVar,icg) = REAL(4*(iVar-1) ,RFREAL) &
563 ! + REAL(4*(iVar-1)+1,RFREAL)*pGrid%cofg(XCOORD,icg) &
564 ! + REAL(4*(iVar-1)+2,RFREAL)*pGrid%cofg(YCOORD,icg) &
565 ! + REAL(4*(iVar-1)+3,RFREAL)*pGrid%cofg(ZCOORD,icg)
566 ! END DO ! iVar
567 ! END DO ! icg
568 ! END DEBUG
569 
570  DO icg = 1,pgrid%nCellsTot
571  r11 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_11)
572  r12 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_12)
573  r22 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_22)
574  r13 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_13)
575  r23 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_23)
576  r33 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_33)
577 
578  c11 = 1.0_rfreal/r11
579  c22 = 1.0_rfreal/r22
580  c33 = 1.0_rfreal/r33
581 
582  c12 = - c11*r12
583  c13 = -(c11*r13 + c12*c22*r23)
584 
585  c23 = - c22*r23
586 
587  cofg(xcoord) = pgrid%cofg(xcoord,icg)
588  cofg(ycoord) = pgrid%cofg(ycoord,icg)
589  cofg(zcoord) = pgrid%cofg(zcoord,icg)
590 
591  DO igrad = ibeggrad,iendgrad
592  grad(xcoord,igrad,icg) = 0.0_rfreal
593  grad(ycoord,igrad,icg) = 0.0_rfreal
594  grad(zcoord,igrad,icg) = 0.0_rfreal
595  END DO ! iVar
596 
597  DO isl = 1,pgrid%c2cs(icg)%nCellMembs
598  icg2 = pgrid%c2cs(icg)%cellMembs(isl)
599 
600  dx = pgrid%cofg(xcoord,icg2) - cofg(xcoord)
601  dy = pgrid%cofg(ycoord,icg2) - cofg(ycoord)
602  dz = pgrid%cofg(zcoord,icg2) - cofg(zcoord)
603 
604  term = 1.0_rfreal/sqrt(dx*dx + dy*dy + dz*dz)
605 
606  dx = term*dx
607  dy = term*dy
608  dz = term*dz
609 
610  term1 = c11*c11*( dx)
611  term2 = c22*c22*( dy + c12*dx)
612  term3 = c33*c33*(dz + c23*dy + c13*dx)
613 
614  wx = term*(term1 + c12*term2 + c13*term3)
615  wy = term*( term2 + c23*term3)
616  wz = term*( term3)
617 
618  igrad = ibeggrad
619 
620  DO ivar = ibegvar,iendvar
621  dvar = var(ivar,icg2) - var(ivar,icg)
622 
623  grad(xcoord,igrad,icg) = grad(xcoord,igrad,icg) + wx*dvar
624  grad(ycoord,igrad,icg) = grad(ycoord,igrad,icg) + wy*dvar
625  grad(zcoord,igrad,icg) = grad(zcoord,igrad,icg) + wz*dvar
626 
627  igrad = igrad + 1
628  END DO ! iVar
629  END DO ! isl
630  END DO ! icg
631 
632 ! DEBUG
633 ! DO iGrad = iBegGrad,iEndGrad
634 ! WRITE(*,*) iGrad,MINVAL(grad(XCOORD,iGrad,1:pGrid%nCellsTot)), &
635 ! MINVAL(grad(YCOORD,iGrad,1:pGrid%nCellsTot)), &
636 ! MINVAL(grad(ZCOORD,iGrad,1:pGrid%nCellsTot)), &
637 ! MAXVAL(grad(XCOORD,iGrad,1:pGrid%nCellsTot)), &
638 ! MAXVAL(grad(YCOORD,iGrad,1:pGrid%nCellsTot)), &
639 ! MAXVAL(grad(ZCOORD,iGrad,1:pGrid%nCellsTot))
640 ! END DO ! iGrad
641 !
642 ! STOP
643 ! END DEBUG
644 
645 ! ******************************************************************************
646 ! End
647 ! ******************************************************************************
648 
649 #ifdef ROCPROF
650  CALL fprofiler_ends("RFLU::ComputeGradCells_3D")
651 #endif
652 
653  CALL deregisterfunction(global)
654 
655  END SUBROUTINE rflu_computegradcells_3d
656 
657 
658 
659 
660 
661 
662 
663 
664 ! ******************************************************************************
665 !
666 ! Purpose: Compute 2D gradients of any vector or scalar at cell centers.
667 !
668 ! Description: None.
669 !
670 ! Input:
671 ! pRegion Pointer to region data
672 ! iBegVar Beginning index of data in var
673 ! iEndVar Ending index of data in var
674 ! iBegGrad Beginning index of data in grad
675 ! iEndGrad Ending index of data in grad
676 ! var Variables of which gradients are to be determined
677 !
678 ! Output:
679 ! grad Gradients of variables at cell centers
680 !
681 ! Notes:
682 ! 1. Use original form of reconstruction, i.e., weights multiplying variable
683 ! differences, because central location is cell, so have variable located
684 ! there. This is in contrast to face gradients and vertex interpolation,
685 ! where need to use modified form of reconstruction.
686 ! 2. If the weighting is changed from inverse-distance to none, then the
687 ! routine RFLU_ComputeStencilMomentsX in RFLU_ModWeights must also be
688 ! adapted.
689 ! 3. Optimized by Adam Moody and Charles Shereda, LLNL.
690 !
691 ! ******************************************************************************
692 
693  SUBROUTINE rflu_computegradcellsfast_2d(pRegion,iBegVar,iEndVar,iBegGrad, &
694  iendgrad,var,grad)
695 
696  IMPLICIT NONE
697 
698 ! ******************************************************************************
699 ! Definitions and declarations
700 ! ******************************************************************************
701 
702 ! ==============================================================================
703 ! Arguments
704 ! ==============================================================================
705 
706  INTEGER :: ibegvar,iendvar,ibeggrad,iendgrad
707  REAL(RFREAL), DIMENSION(:,:), POINTER :: var
708  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
709  TYPE(t_region), POINTER :: pregion
710 
711 ! ==============================================================================
712 ! Locals
713 ! ==============================================================================
714 
715  INTEGER :: endisl,errorflag,icg,icg2,icg3,icg4,igrad,isl,ivar
716  REAL(RFREAL) :: c11,c1111,c12,c22,c2222,dvar,dvar3,dvar4,dx,dx_icg,dy, &
717  dy_icg,r11,r12,r22,term,term1,term2,wtx,wty
718  TYPE(t_global), POINTER :: global
719  TYPE(t_grid), POINTER :: pgrid
720 
721 ! *****************************************************************************
722 ! Start
723 ! *****************************************************************************
724 
725  global => pregion%global
726 
727  CALL registerfunction(global,'RFLU_ComputeGradCellsFast_2D',&
728  'RFLU_ModDifferentiationCells.F90' )
729 
730 #ifdef ROCPROF
731  CALL fprofiler_begins("RFLU::ComputeGradCellsFast_2D")
732 #endif
733 
734  IF ( (iendvar - ibegvar) /= (iendgrad - ibeggrad) ) THEN
735  CALL errorstop(global,err_grad_mismatch,__line__)
736  END IF ! iEndVar
737 
738 ! ******************************************************************************
739 ! Set pointers and variables
740 ! ******************************************************************************
741 
742  pgrid => pregion%grid
743 
744 ! ******************************************************************************
745 ! Loop over cells and compute gradients
746 ! ******************************************************************************
747 
748 ! DEBUG
749 ! DO icg = 1,pGrid%nCellsTot
750 ! DO iVar = iBegVar,iEndVar
751 ! var(iVar,icg) = REAL(4*(iVar-1) ,RFREAL) &
752 ! + REAL(4*(iVar-1)+1,RFREAL)*pGrid%cofg(XCOORD,icg) &
753 ! + REAL(4*(iVar-1)+2,RFREAL)*pGrid%cofg(YCOORD,icg)
754 ! END DO ! iVar
755 ! END DO ! icg
756 ! END DEBUG
757 
758  DO icg = 1,pgrid%nCellsTot
759  r11 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_11)
760  r12 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_12)
761  r22 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_22)
762 
763  c11 = 1.0_rfreal/r11
764  c22 = 1.0_rfreal/r22
765 
766  c12 = -c11*r12
767 
768  endisl = pgrid%c2cs(icg)%nCellMembs
769 
770  DO igrad = ibeggrad,iendgrad
771  grad(xcoord,igrad,icg) = 0.0_rfreal
772  grad(ycoord,igrad,icg) = 0.0_rfreal
773  grad(zcoord,igrad,icg) = 0.0_rfreal ! Set z-component to zero
774  END DO ! iGrad
775 
776  icg3 = pgrid%c2cs(icg)%cellMembs(1)
777  icg4 = pgrid%c2cs(icg)%cellMembs(2)
778 
779  dx_icg = pgrid%cofg(xcoord,icg)
780  dy_icg = pgrid%cofg(ycoord,icg)
781 
782  dvar3 = var(ibegvar,icg3) - var(ibegvar,icg)
783  dvar4 = var(ibegvar,icg4) - var(ibegvar,icg)
784 
785  c1111 = c11*c11
786  c2222 = c22*c22
787 
788  DO isl = 1,endisl-2
789  icg2 = icg3
790  icg3 = icg4
791  icg4 = pgrid%c2cs(icg)%cellMembs(isl+2)
792 
793  dvar = dvar3
794  dvar3 = dvar4
795  dvar4 = var(ibegvar,icg4) - var(ibegvar,icg)
796 
797  dx = pgrid%cofg(xcoord,icg2) - dx_icg
798  dy = pgrid%cofg(ycoord,icg2) - dy_icg
799 
800  term = 1.0_rfreal/sqrt(dx*dx + dy*dy)
801 
802  dx = term*dx
803  dy = term*dy
804 
805  term1 = c1111*( dx)
806  term2 = c2222*(dy + c12*dx)
807 
808  wtx = term*(term1 + c12*term2)
809  wty = term*( term2)
810 
811  igrad = ibeggrad
812 
813  grad(xcoord,igrad,icg) = grad(xcoord,igrad,icg) + wtx*dvar
814  grad(ycoord,igrad,icg) = grad(ycoord,igrad,icg) + wty*dvar
815 
816  igrad = igrad + 1
817 
818  DO ivar = ibegvar+1,iendvar
819  dvar = var(ivar,icg2) - var(ivar,icg)
820 
821  grad(xcoord,igrad,icg) = grad(xcoord,igrad,icg) + wtx*dvar
822  grad(ycoord,igrad,icg) = grad(ycoord,igrad,icg) + wty*dvar
823 
824  igrad = igrad + 1
825  END DO ! iVar
826  END DO ! isl
827 
828  DO isl = endisl-1,endisl
829  icg2 = icg3
830  icg3 = icg4
831 
832  dvar = dvar3
833  dvar3 = dvar4
834 
835  dx = pgrid%cofg(xcoord,icg2) - dx_icg
836  dy = pgrid%cofg(ycoord,icg2) - dy_icg
837 
838  term = 1.0_rfreal/sqrt(dx*dx + dy*dy)
839 
840  dx = term*dx
841  dy = term*dy
842 
843  term1 = c1111*( dx)
844  term2 = c2222*(dy + c12*dx)
845 
846  wtx = term*(term1 + c12*term2)
847  wty = term*( term2)
848 
849  igrad = ibeggrad
850 
851  grad(xcoord,igrad,icg) = grad(xcoord,igrad,icg) + wtx*dvar
852  grad(ycoord,igrad,icg) = grad(ycoord,igrad,icg) + wty*dvar
853 
854  igrad = igrad + 1
855 
856  DO ivar = ibegvar+1,iendvar
857  dvar = var(ivar,icg2) - var(ivar,icg)
858 
859  grad(xcoord,igrad,icg) = grad(xcoord,igrad,icg) + wtx*dvar
860  grad(ycoord,igrad,icg) = grad(ycoord,igrad,icg) + wty*dvar
861 
862  igrad = igrad + 1
863  END DO ! iVar
864  END DO ! isl
865  END DO ! icg
866 
867 ! DEBUG
868 ! DO iGrad = iBegGrad,iEndGrad
869 ! WRITE(*,*) iGrad,MINVAL(grad(XCOORD,iGrad,1:pGrid%nCellsTot)), &
870 ! MINVAL(grad(YCOORD,iGrad,1:pGrid%nCellsTot)), &
871 ! MAXVAL(grad(XCOORD,iGrad,1:pGrid%nCellsTot)), &
872 ! MAXVAL(grad(YCOORD,iGrad,1:pGrid%nCellsTot))
873 ! END DO ! iGrad
874 !
875 ! STOP
876 ! END DEBUG
877 
878 ! ******************************************************************************
879 ! End
880 ! ******************************************************************************
881 
882 #ifdef ROCPROF
883  CALL fprofiler_ends("RFLU::ComputeGradCellsFast_2D")
884 #endif
885 
886  CALL deregisterfunction(global)
887 
888  END SUBROUTINE rflu_computegradcellsfast_2d
889 
890 
891 
892 
893 
894 
895 
896 ! ******************************************************************************
897 !
898 ! Purpose: Compute 3D gradients of any vector or scalar at cell centers.
899 !
900 ! Description: None.
901 !
902 ! Input:
903 ! pRegion Pointer to region data
904 ! iBegVar Beginning index of data in var
905 ! iEndVar Ending index of data in var
906 ! iBegGrad Beginning index of data in grad
907 ! iEndGrad Ending index of data in grad
908 ! var Variables of which gradients are to be determined
909 !
910 ! Output:
911 ! grad Gradients of variables at cell centers
912 !
913 ! Notes:
914 ! 1. Use original form of reconstruction, i.e., weights multiplying variable
915 ! differences, because central location is cell, so have variable located
916 ! there. This is in contrast to face gradients and vertex interpolation,
917 ! where need to use modified form of reconstruction.
918 ! 2. If the weighting is changed from inverse-distance to none, then the
919 ! routine RFLU_ComputeStencilMomentsX in RFLU_ModWeights must also be
920 ! adapted.
921 ! 3. Optimized by Adam Moody and Charles Shereda, LLNL.
922 !
923 ! ******************************************************************************
924 
925  SUBROUTINE rflu_computegradcellsfast_3d(pRegion,iBegVar,iEndVar,iBegGrad, &
926  iendgrad,var,grad)
927 
928  IMPLICIT NONE
929 
930 ! ******************************************************************************
931 ! Definitions and declarations
932 ! ******************************************************************************
933 
934 ! ==============================================================================
935 ! Arguments
936 ! ==============================================================================
937 
938  INTEGER :: ibegvar,iendvar,ibeggrad,iendgrad
939  REAL(RFREAL), DIMENSION(:,:), POINTER :: var
940  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
941  TYPE(t_region), POINTER :: pregion
942 
943 ! ==============================================================================
944 ! Locals
945 ! ==============================================================================
946 
947  INTEGER :: endisl,errorflag,icg,icg2,icg3,icg4,igrad,isl,ivar
948  REAL(RFREAL) :: c11,c1111,c12,c13,c22,c2222,c23,c33,c3333,dvar,dvar3, &
949  dvar4,dx,dx_icg,dy,dy_icg,dz,dz_icg,r11,r12,r13,r22, &
950  r23,r33,term,term1,term2,term3,wtx,wty,wtz
951  TYPE(t_global), POINTER :: global
952  TYPE(t_grid), POINTER :: pgrid
953 
954 ! *****************************************************************************
955 ! Start
956 ! *****************************************************************************
957 
958  global => pregion%global
959 
960  CALL registerfunction(global,'RFLU_ComputeGradCellsFast_3D',&
961  'RFLU_ModDifferentiationCells.F90' )
962 
963 #ifdef ROCPROF
964  CALL fprofiler_begins("RFLU::ComputeGradCellsFast_3D")
965 #endif
966 
967  IF ( (iendvar - ibegvar) /= (iendgrad - ibeggrad) ) THEN
968  CALL errorstop(global,err_grad_mismatch,__line__)
969  END IF ! iEndVar
970 
971 ! ******************************************************************************
972 ! Set pointers and variables
973 ! ******************************************************************************
974 
975  pgrid => pregion%grid
976 
977 ! ******************************************************************************
978 ! Loop over cells and compute gradients
979 ! ******************************************************************************
980 
981 ! DEBUG
982 ! DO icg = 1,pGrid%nCellsTot
983 ! DO iVar = iBegVar,iEndVar
984 ! var(iVar,icg) = REAL(4*(iVar-1) ,RFREAL) &
985 ! + REAL(4*(iVar-1)+1,RFREAL)*pGrid%cofg(XCOORD,icg) &
986 ! + REAL(4*(iVar-1)+2,RFREAL)*pGrid%cofg(YCOORD,icg) &
987 ! + REAL(4*(iVar-1)+3,RFREAL)*pGrid%cofg(ZCOORD,icg)
988 ! END DO ! iVar
989 ! END DO ! icg
990 ! END DEBUG
991 
992  DO icg = 1,pgrid%nCellsTot
993  r11 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_11)
994  r12 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_12)
995  r22 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_22)
996  r13 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_13)
997  r23 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_23)
998  r33 = pgrid%c2cs(icg)%xyzMoms(xyz_mom_33)
999 
1000  c11 = 1.0_rfreal/r11
1001  c22 = 1.0_rfreal/r22
1002  c33 = 1.0_rfreal/r33
1003 
1004  c12 = - c11*r12
1005  c13 = -(c11*r13 + c12*c22*r23)
1006 
1007  c23 = - c22*r23
1008 
1009  endisl = pgrid%c2cs(icg)%nCellMembs
1010 
1011  DO igrad = ibeggrad,iendgrad
1012  grad(xcoord,igrad,icg) = 0.0_rfreal
1013  grad(ycoord,igrad,icg) = 0.0_rfreal
1014  grad(zcoord,igrad,icg) = 0.0_rfreal
1015  END DO ! iGrad
1016 
1017  icg3 = pgrid%c2cs(icg)%cellMembs(1)
1018  icg4 = pgrid%c2cs(icg)%cellMembs(2)
1019 
1020  dx_icg = pgrid%cofg(xcoord,icg)
1021  dy_icg = pgrid%cofg(ycoord,icg)
1022  dz_icg = pgrid%cofg(zcoord,icg)
1023 
1024  dvar3 = var(ibegvar,icg3) - var(ibegvar,icg)
1025  dvar4 = var(ibegvar,icg4) - var(ibegvar,icg)
1026 
1027  c1111 = c11*c11
1028  c2222 = c22*c22
1029  c3333 = c33*c33
1030 
1031  DO isl = 1,endisl-2
1032  icg2 = icg3
1033  icg3 = icg4
1034  icg4 = pgrid%c2cs(icg)%cellMembs(isl+2)
1035 
1036  dvar = dvar3
1037  dvar3 = dvar4
1038  dvar4 = var(ibegvar,icg4) - var(ibegvar,icg)
1039 
1040  dx = pgrid%cofg(xcoord,icg2) - dx_icg
1041  dy = pgrid%cofg(ycoord,icg2) - dy_icg
1042  dz = pgrid%cofg(zcoord,icg2) - dz_icg
1043 
1044  term = 1.0_rfreal/sqrt(dx*dx + dy*dy + dz*dz)
1045 
1046  dx = term*dx
1047  dy = term*dy
1048  dz = term*dz
1049 
1050  term1 = c1111*( dx)
1051  term2 = c2222*( dy + c12*dx)
1052  term3 = c3333*(dz + c23*dy + c13*dx)
1053 
1054  wtx = term*(term1 + c12*term2 + c13*term3)
1055  wty = term*( term2 + c23*term3)
1056  wtz = term*( term3)
1057 
1058  igrad = ibeggrad
1059 
1060  grad(xcoord,igrad,icg) = grad(xcoord,igrad,icg) + wtx*dvar
1061  grad(ycoord,igrad,icg) = grad(ycoord,igrad,icg) + wty*dvar
1062  grad(zcoord,igrad,icg) = grad(zcoord,igrad,icg) + wtz*dvar
1063 
1064  igrad = igrad + 1
1065 
1066  DO ivar = ibegvar+1,iendvar
1067  dvar = var(ivar,icg2) - var(ivar,icg)
1068 
1069  grad(xcoord,igrad,icg) = grad(xcoord,igrad,icg) + wtx*dvar
1070  grad(ycoord,igrad,icg) = grad(ycoord,igrad,icg) + wty*dvar
1071  grad(zcoord,igrad,icg) = grad(zcoord,igrad,icg) + wtz*dvar
1072 
1073  igrad = igrad + 1
1074  END DO ! iVar
1075  END DO ! isl
1076 
1077  DO isl = endisl-1,endisl
1078  icg2 = icg3
1079  icg3 = icg4
1080 
1081  dvar = dvar3
1082  dvar3 = dvar4
1083 
1084  dx = pgrid%cofg(xcoord,icg2) - dx_icg
1085  dy = pgrid%cofg(ycoord,icg2) - dy_icg
1086  dz = pgrid%cofg(zcoord,icg2) - dz_icg
1087 
1088  term = 1.0_rfreal/sqrt(dx*dx + dy*dy + dz*dz)
1089 
1090  dx = term*dx
1091  dy = term*dy
1092  dz = term*dz
1093 
1094  term1 = c1111*( dx)
1095  term2 = c2222*( dy + c12*dx)
1096  term3 = c3333*(dz + c23*dy + c13*dx)
1097 
1098  wtx = term*(term1 + c12*term2 + c13*term3)
1099  wty = term*( term2 + c23*term3)
1100  wtz = term*( term3)
1101 
1102  igrad = ibeggrad
1103 
1104  grad(xcoord,igrad,icg) = grad(xcoord,igrad,icg) + wtx*dvar
1105  grad(ycoord,igrad,icg) = grad(ycoord,igrad,icg) + wty*dvar
1106  grad(zcoord,igrad,icg) = grad(zcoord,igrad,icg) + wtz*dvar
1107 
1108  igrad = igrad + 1
1109 
1110  DO ivar = ibegvar+1,iendvar
1111  dvar = var(ivar,icg2) - var(ivar,icg)
1112 
1113  grad(xcoord,igrad,icg) = grad(xcoord,igrad,icg) + wtx*dvar
1114  grad(ycoord,igrad,icg) = grad(ycoord,igrad,icg) + wty*dvar
1115  grad(zcoord,igrad,icg) = grad(zcoord,igrad,icg) + wtz*dvar
1116 
1117  igrad = igrad + 1
1118  END DO ! iVar
1119  END DO ! isl
1120  END DO ! icg
1121 
1122 ! DEBUG
1123 ! DO iGrad = iBegGrad,iEndGrad
1124 ! WRITE(*,*) iGrad,MINVAL(grad(XCOORD,iGrad,1:pGrid%nCellsTot)), &
1125 ! MINVAL(grad(YCOORD,iGrad,1:pGrid%nCellsTot)), &
1126 ! MINVAL(grad(ZCOORD,iGrad,1:pGrid%nCellsTot)), &
1127 ! MAXVAL(grad(XCOORD,iGrad,1:pGrid%nCellsTot)), &
1128 ! MAXVAL(grad(YCOORD,iGrad,1:pGrid%nCellsTot)), &
1129 ! MAXVAL(grad(ZCOORD,iGrad,1:pGrid%nCellsTot))
1130 ! END DO ! iGrad
1131 !
1132 ! STOP
1133 ! END DEBUG
1134 
1135 ! ******************************************************************************
1136 ! End
1137 ! ******************************************************************************
1138 
1139 #ifdef ROCPROF
1140  CALL fprofiler_ends("RFLU::ComputeGradCellsFast_3D")
1141 #endif
1142 
1143  CALL deregisterfunction(global)
1144 
1145  END SUBROUTINE rflu_computegradcellsfast_3d
1146 
1147 
1148 
1149 
1150 
1151 
1152 
1153 ! ******************************************************************************
1154 !
1155 ! Purpose: Compute constrained gradients of any vector or scalar at cell
1156 ! centers.
1157 !
1158 ! Description: None.
1159 !
1160 ! Input:
1161 ! pRegion Pointer to region data
1162 ! iBegVar Beginning index of data in var
1163 ! iEndVar Ending index of data in var
1164 ! iBegGrad Beginning index of data in grad
1165 ! iEndGrad Ending index of data in grad
1166 ! varInfo Variable information
1167 ! var Variables of which gradients are to be determined
1168 !
1169 ! Output:
1170 ! grad Gradients of variables at cell centers
1171 !
1172 ! Notes:
1173 ! 1. If the weighting is changed from inverse-distance to none, then the
1174 ! routine RFLU_ComputeStencilMomentsX in RFLU_ModWeights must also be
1175 ! adapted.
1176 ! 2. Restricted to linear reconstruction for now.
1177 ! 3. Restricted to Dirichlet constraints for now.
1178 !
1179 ! ******************************************************************************
1180 
1181  SUBROUTINE rflu_computegradcellsconstr(pRegion,iBegVar,iEndVar,iBegGrad, &
1182  iendgrad,varinfo,var,grad)
1183 
1185 
1186  IMPLICIT NONE
1187 
1188 ! ******************************************************************************
1189 ! Definitions and declarations
1190 ! ******************************************************************************
1191 
1192 ! ==============================================================================
1193 ! Arguments
1194 ! ==============================================================================
1195 
1196  INTEGER, INTENT(IN) :: ibegvar,iendvar,ibeggrad,iendgrad
1197  INTEGER, INTENT(IN) :: varinfo(ibegvar:iendvar)
1198  REAL(RFREAL), DIMENSION(:,:), POINTER :: var
1199  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
1200  TYPE(t_region), POINTER :: pregion
1201 
1202 ! ==============================================================================
1203 ! Locals
1204 ! ==============================================================================
1205 
1206  INTEGER :: errorflag,icg,icg2,icl,icol,ifg,ifl,igrad,ipatch,isl,irow,ivar, &
1207  ncols,nconstr,nrows,scount
1208  INTEGER, DIMENSION(:), ALLOCATABLE :: constrtype
1209  REAL(RFREAL) :: cwt,dvar,dx,dy,dz,term,varc,gx,gy,gz
1210  REAL(RFREAL) :: cofg(xcoord:zcoord)
1211  REAL(RFREAL) :: colmax(4)
1212  REAL(RFREAL), DIMENSION(:,:), ALLOCATABLE :: a,ainv
1213  TYPE(t_global), POINTER :: global
1214  TYPE(t_grid), POINTER :: pgrid
1215  TYPE(t_patch), POINTER :: ppatch
1216 
1217 ! ******************************************************************************
1218 ! Start
1219 ! ******************************************************************************
1220 
1221  global => pregion%global
1222 
1223  CALL registerfunction(global,'RFLU_ComputeGradCellsConstr',&
1224  'RFLU_ModDifferentiationCells.F90')
1225 
1226 #ifdef ROCPROF
1227  CALL fprofiler_begins("RFLU::ComputeGradCellsConstr")
1228 #endif
1229 
1230  IF ( (iendvar - ibegvar) /= (iendgrad - ibeggrad) ) THEN
1231  CALL errorstop(global,err_grad_mismatch,__line__)
1232  END IF ! iEndVar
1233 
1234 ! ******************************************************************************
1235 ! Set pointers and variables
1236 ! ******************************************************************************
1237 
1238  pgrid => pregion%grid
1239 
1240  cwt = pregion%mixtInput%cReconstCellsWeight
1241 
1242 ! ******************************************************************************
1243 ! Loop over cells and compute constrained gradients
1244 ! ******************************************************************************
1245 
1246  DO icl = 1,pgrid%nCellsConstr
1247  icg = pgrid%icgConstr(icl)
1248 
1249 ! ==============================================================================
1250 ! Initialize gradients, set cofg, and allocate memory for constraint type
1251 ! ==============================================================================
1252 
1253  DO igrad = ibeggrad,iendgrad
1254  grad(xcoord,igrad,icg) = 0.0_rfreal
1255  grad(ycoord,igrad,icg) = 0.0_rfreal
1256  grad(zcoord,igrad,icg) = 0.0_rfreal
1257  END DO ! iGrad
1258 
1259  cofg(xcoord) = pgrid%cofg(xcoord,icg)
1260  cofg(ycoord) = pgrid%cofg(ycoord,icg)
1261  cofg(zcoord) = pgrid%cofg(zcoord,icg)
1262 
1263  ALLOCATE(constrtype(pgrid%c2cs(icg)%nBFaceMembs),stat=errorflag)
1264  global%error = errorflag
1265  IF ( global%error /= err_none ) THEN
1266  CALL errorstop(global,err_allocate,__line__,'constrType')
1267  END IF ! global%error
1268 
1269 ! ==============================================================================
1270 ! Compute gradients
1271 ! ==============================================================================
1272 
1273  igrad = ibeggrad
1274 
1275  DO ivar = ibegvar,iendvar
1276 
1277 ! ------------------------------------------------------------------------------
1278 ! Determine number of constraints
1279 ! ------------------------------------------------------------------------------
1280 
1281  nconstr = 0
1282 
1283  DO isl = 1,pgrid%c2cs(icg)%nBFaceMembs
1284  ipatch = pgrid%c2cs(icg)%bFaceMembs(1,isl)
1285  ifl = pgrid%c2cs(icg)%bFaceMembs(2,isl)
1286 
1287  ppatch => pregion%patches(ipatch)
1288 
1289  constrtype(isl) = rflu_getconstrtype(pregion,ppatch,varinfo(ivar),ifl)
1290 
1291  IF ( constrtype(isl) == constr_type_dirichlet ) THEN
1292  nconstr = nconstr + 1
1293  ELSE
1294  constrtype(isl) = constr_type_none
1295  END IF ! constrType
1296  END DO ! isl
1297 
1298 ! ------------------------------------------------------------------------------
1299 ! Gradients constrained by Dirichlet boundary conditions. Treated as
1300 ! soft constraints. NOTE do not need to treat case of unconstrained
1301 ! gradients here because they should already have been computed. If
1302 ! they have not been computed at this stage, they will be set to zero
1303 ! by initialization above.
1304 ! ------------------------------------------------------------------------------
1305 
1306  IF ( nconstr > 0 ) THEN
1307 
1308 ! ------- Allocate temporary memory --------------------------------------------
1309 
1310  nrows = pgrid%c2cs(icg)%nCellMembs + nconstr
1311  ncols = pregion%mixtInput%dimens
1312 
1313  ALLOCATE(a(nrows,ncols),stat=errorflag)
1314  global%error = errorflag
1315  IF ( global%error /= err_none ) THEN
1316  CALL errorstop(global,err_allocate,__line__,'a')
1317  END IF ! global%error
1318 
1319  ALLOCATE(ainv(ncols,nrows),stat=errorflag)
1320  global%error = errorflag
1321  IF ( global%error /= err_none ) THEN
1322  CALL errorstop(global,err_allocate,__line__,'aInv')
1323  END IF ! global%error
1324 
1325 ! ------- Define left-hand side matrix -----------------------------------------
1326 
1327  SELECT CASE ( pregion%mixtInput%dimens )
1328 
1329 ! --------- Two dimensions
1330 
1331  CASE ( 2 )
1332  DO isl = 1,pgrid%c2cs(icg)%nCellMembs
1333  icg2 = pgrid%c2cs(icg)%cellMembs(isl)
1334 
1335  dx = pgrid%cofg(xcoord,icg2) - cofg(xcoord)
1336  dy = pgrid%cofg(ycoord,icg2) - cofg(ycoord)
1337 
1338  term = 1.0_rfreal/sqrt(dx*dx + dy*dy)
1339 
1340  a(isl,1) = term*dx
1341  a(isl,2) = term*dy
1342  END DO ! isl
1343 
1344  irow = pgrid%c2cs(icg)%nCellMembs
1345 
1346  DO isl = 1,pgrid%c2cs(icg)%nBFaceMembs
1347  IF ( constrtype(isl) == constr_type_dirichlet ) THEN
1348  ipatch = pgrid%c2cs(icg)%bFaceMembs(1,isl)
1349  ifl = pgrid%c2cs(icg)%bFaceMembs(2,isl)
1350 
1351  ppatch => pregion%patches(ipatch)
1352 
1353  dx = ppatch%fc(xcoord,ifl) - cofg(xcoord)
1354  dy = ppatch%fc(ycoord,ifl) - cofg(ycoord)
1355 
1356  term = cwt/sqrt(dx*dx + dy*dy)
1357 
1358  irow = irow + 1
1359 
1360  a(irow,1) = term*dx
1361  a(irow,2) = term*dy
1362  END IF ! constrType
1363  END DO ! isl
1364 
1365 ! --------- Three dimensions
1366 
1367  CASE ( 3 )
1368  DO isl = 1,pgrid%c2cs(icg)%nCellMembs
1369  icg2 = pgrid%c2cs(icg)%cellMembs(isl)
1370 
1371  dx = pgrid%cofg(xcoord,icg2) - cofg(xcoord)
1372  dy = pgrid%cofg(ycoord,icg2) - cofg(ycoord)
1373  dz = pgrid%cofg(zcoord,icg2) - cofg(zcoord)
1374 
1375  term = 1.0_rfreal/sqrt(dx*dx + dy*dy + dz*dz)
1376 
1377  a(isl,1) = term*dx
1378  a(isl,2) = term*dy
1379  a(isl,3) = term*dz
1380  END DO ! isl
1381 
1382  irow = pgrid%c2cs(icg)%nCellMembs
1383 
1384  DO isl = 1,pgrid%c2cs(icg)%nBFaceMembs
1385  IF ( constrtype(isl) == constr_type_dirichlet ) THEN
1386  ipatch = pgrid%c2cs(icg)%bFaceMembs(1,isl)
1387  ifl = pgrid%c2cs(icg)%bFaceMembs(2,isl)
1388 
1389  ppatch => pregion%patches(ipatch)
1390 
1391  dx = ppatch%fc(xcoord,ifl) - cofg(xcoord)
1392  dy = ppatch%fc(ycoord,ifl) - cofg(ycoord)
1393  dz = ppatch%fc(zcoord,ifl) - cofg(zcoord)
1394 
1395  term = cwt/sqrt(dx*dx + dy*dy + dz*dz)
1396 
1397  irow = irow + 1
1398 
1399  a(irow,1) = term*dx
1400  a(irow,2) = term*dy
1401  a(irow,3) = term*dz
1402  END IF ! constrType
1403  END DO ! isl
1404  CASE default
1405  CALL errorstop(global,err_reached_default,__line__)
1406  END SELECT ! pRegion%mixtInput%dimens
1407 
1408 ! ------- Compute constrained gradient weights ---------------------------------
1409 
1410  DO icol = 1,ncols
1411  colmax(icol) = -huge(1.0_rfreal)
1412 
1413  DO irow = 1,nrows
1414  colmax(icol) = max(colmax(icol),abs(a(irow,icol)))
1415  END DO ! iRow
1416 
1417  DO irow = 1,nrows
1418  a(irow,icol) = a(irow,icol)/colmax(icol)
1419  END DO ! iRow
1420  END DO ! iCol
1421 
1422  CALL rflu_invertmatrixsvd(global,nrows,ncols,a,ainv,scount)
1423 
1424  DO icol = 1,ncols
1425  DO irow = 1,nrows
1426  ainv(icol,irow) = ainv(icol,irow)/colmax(icol)
1427  END DO ! iRow
1428  END DO ! iCol
1429 
1430 ! TEMPORARY
1431  IF ( scount /= 0 ) THEN
1432  WRITE(*,*) 'ERROR - Singular matrix in RFLU_ComputeGradCellsConstr!'
1433  stop
1434  END IF ! sCount
1435 ! END TEMPORARY
1436 
1437 ! ------- Compute constrained gradients ----------------------------------------
1438 
1439  SELECT CASE ( pregion%mixtInput%dimens )
1440 
1441 ! --------- Two dimensions
1442 
1443  CASE ( 2 )
1444  gx = 0.0_rfreal
1445  gy = 0.0_rfreal
1446 
1447  DO isl = 1,pgrid%c2cs(icg)%nCellMembs
1448  icg2 = pgrid%c2cs(icg)%cellMembs(isl)
1449 
1450  dx = pgrid%cofg(xcoord,icg2) - cofg(xcoord)
1451  dy = pgrid%cofg(ycoord,icg2) - cofg(ycoord)
1452 
1453  term = 1.0_rfreal/sqrt(dx*dx + dy*dy)
1454 
1455  dvar = var(ivar,icg2) - var(ivar,icg)
1456 
1457  gx = gx + term*ainv(1,isl)*dvar
1458  gy = gy + term*ainv(2,isl)*dvar
1459  END DO ! isl
1460 
1461  irow = pgrid%c2cs(icg)%nCellMembs
1462 
1463  DO isl = 1,pgrid%c2cs(icg)%nBFaceMembs
1464  IF ( constrtype(isl) == constr_type_dirichlet ) THEN
1465  ipatch = pgrid%c2cs(icg)%bFaceMembs(1,isl)
1466  ifl = pgrid%c2cs(icg)%bFaceMembs(2,isl)
1467 
1468  ppatch => pregion%patches(ipatch)
1469 
1470  dx = ppatch%fc(xcoord,ifl) - cofg(xcoord)
1471  dy = ppatch%fc(ycoord,ifl) - cofg(ycoord)
1472 
1473  term = cwt/sqrt(dx*dx + dy*dy)
1474 
1475  varc = rflu_getconstrvalue(pregion,ppatch,varinfo(ivar),ifl)
1476  dvar = varc - var(ivar,icg)
1477 
1478  irow = irow + 1
1479 
1480  gx = gx + term*ainv(1,irow)*dvar
1481  gy = gy + term*ainv(2,irow)*dvar
1482  END IF ! constrType
1483  END DO ! isl
1484 
1485  grad(xcoord,igrad,icg) = gx
1486  grad(ycoord,igrad,icg) = gy
1487  grad(zcoord,igrad,icg) = 0.0_rfreal
1488 
1489 ! --------- Three dimensions
1490 
1491  CASE ( 3 )
1492  gx = 0.0_rfreal
1493  gy = 0.0_rfreal
1494  gz = 0.0_rfreal
1495 
1496  DO isl = 1,pgrid%c2cs(icg)%nCellMembs
1497  icg2 = pgrid%c2cs(icg)%cellMembs(isl)
1498 
1499  dx = pgrid%cofg(xcoord,icg2) - cofg(xcoord)
1500  dy = pgrid%cofg(ycoord,icg2) - cofg(ycoord)
1501  dz = pgrid%cofg(zcoord,icg2) - cofg(zcoord)
1502 
1503  term = 1.0_rfreal/sqrt(dx*dx + dy*dy + dz*dz)
1504 
1505  dvar = var(ivar,icg2) - var(ivar,icg)
1506 
1507  gx = gx + term*ainv(1,isl)*dvar
1508  gy = gy + term*ainv(2,isl)*dvar
1509  gz = gz + term*ainv(3,isl)*dvar
1510  END DO ! isl
1511 
1512  irow = pgrid%c2cs(icg)%nCellMembs
1513 
1514  DO isl = 1,pgrid%c2cs(icg)%nBFaceMembs
1515  IF ( constrtype(isl) == constr_type_dirichlet ) THEN
1516  ipatch = pgrid%c2cs(icg)%bFaceMembs(1,isl)
1517  ifl = pgrid%c2cs(icg)%bFaceMembs(2,isl)
1518 
1519  ppatch => pregion%patches(ipatch)
1520 
1521  dx = ppatch%fc(xcoord,ifl) - cofg(xcoord)
1522  dy = ppatch%fc(ycoord,ifl) - cofg(ycoord)
1523  dz = ppatch%fc(zcoord,ifl) - cofg(zcoord)
1524 
1525  term = cwt/sqrt(dx*dx + dy*dy + dz*dz)
1526 
1527  varc = rflu_getconstrvalue(pregion,ppatch,varinfo(ivar),ifl)
1528  dvar = varc - var(ivar,icg)
1529 
1530  irow = irow + 1
1531 
1532  gx = gx + term*ainv(1,irow)*dvar
1533  gy = gy + term*ainv(2,irow)*dvar
1534  gz = gz + term*ainv(3,irow)*dvar
1535  END IF ! constrType
1536  END DO ! isl
1537 
1538  grad(xcoord,igrad,icg) = gx
1539  grad(ycoord,igrad,icg) = gy
1540  grad(zcoord,igrad,icg) = gz
1541  CASE default
1542  CALL errorstop(global,err_reached_default,__line__)
1543  END SELECT ! pRegion%mixtInput%dimens
1544 
1545 ! ------- Deallocate temporary memory ------------------------------------------
1546 
1547  DEALLOCATE(a,stat=errorflag)
1548  global%error = errorflag
1549  IF ( global%error /= err_none ) THEN
1550  CALL errorstop(global,err_deallocate,__line__,'a')
1551  END IF ! global%error
1552 
1553  DEALLOCATE(ainv,stat=errorflag)
1554  global%error = errorflag
1555  IF ( global%error /= err_none ) THEN
1556  CALL errorstop(global,err_deallocate,__line__,'aInv')
1557  END IF ! global%error
1558  END IF ! gradType
1559 
1560  igrad = igrad + 1
1561  END DO ! iVar
1562 
1563  DEALLOCATE(constrtype,stat=errorflag)
1564  global%error = errorflag
1565  IF ( global%error /= err_none ) THEN
1566  CALL errorstop(global,err_deallocate,__line__,'constrType')
1567  END IF ! global%error
1568  END DO ! icl
1569 
1570 ! ******************************************************************************
1571 ! End
1572 ! ******************************************************************************
1573 
1574 #ifdef ROCPROF
1575  CALL fprofiler_ends("RFLU::ComputeGradCellsConstr")
1576 #endif
1577 
1578  CALL deregisterfunction(global)
1579 
1580  END SUBROUTINE rflu_computegradcellsconstr
1581 
1582 
1583 
1584 
1585 
1586 
1587 
1588 
1589 
1590 ! ******************************************************************************
1591 !
1592 ! Purpose: Compute gradients of any vector or scalar at cell centers.
1593 !
1594 ! Description: None.
1595 !
1596 ! Input:
1597 ! pRegion Pointer to region data
1598 ! iBegVar Beginning index of data in var
1599 ! iEndVar Ending index of data in var
1600 ! iBegGrad Beginning index of data in grad
1601 ! iEndGrad Ending index of data in grad
1602 ! var Variables of which gradients are to be determined
1603 ! varInfo Variable information
1604 !
1605 ! Output:
1606 ! grad Gradients of variables at cell centers
1607 !
1608 ! Notes: None.
1609 !
1610 ! ******************************************************************************
1611 
1612  SUBROUTINE rflu_computegradcellswrapper(pRegion,iBegVar,iEndVar,iBegGrad, &
1613  iendgrad,varinfo,var,grad)
1614 
1615  IMPLICIT NONE
1616 
1617 ! ******************************************************************************
1618 ! Definitions and declarations
1619 ! ******************************************************************************
1620 
1621 ! ==============================================================================
1622 ! Arguments
1623 ! ==============================================================================
1624 
1625  INTEGER, INTENT(IN) :: ibegvar,iendvar,ibeggrad,iendgrad
1626  INTEGER, INTENT(IN) :: varinfo(ibegvar:iendvar)
1627  REAL(RFREAL), DIMENSION(:,:), POINTER :: var
1628  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
1629  TYPE(t_region), POINTER :: pregion
1630 
1631 ! ==============================================================================
1632 ! Locals
1633 ! ==============================================================================
1634 
1635  TYPE(t_global), POINTER :: global
1636  TYPE(t_grid), POINTER :: pgrid
1637 
1638 ! *****************************************************************************
1639 ! Start
1640 ! *****************************************************************************
1641 
1642  global => pregion%global
1643 
1644  CALL registerfunction(global,'RFLU_ComputeGradCellsWrapper',&
1645  'RFLU_ModDifferentiationCells.F90' )
1646 
1647 ! ******************************************************************************
1648 ! Call gradient routines
1649 ! ******************************************************************************
1650 
1651  SELECT CASE ( pregion%mixtInput%stencilDimensCells )
1652  CASE ( 1 )
1653  CALL rflu_computegradcells_1d(pregion,ibegvar,iendvar,ibeggrad, &
1654  iendgrad,var,grad)
1655  CASE ( 2 )
1656 ! CALL RFLU_ComputeGradCells_2D(pRegion,iBegVar,iEndVar,iBegGrad, &
1657 ! iEndGrad,var,grad)
1658  CALL rflu_computegradcellsfast_2d(pregion,ibegvar,iendvar,ibeggrad, &
1659  iendgrad,var,grad)
1660 
1661  IF ( pregion%grid%nCellsConstr > 0 ) THEN
1662  CALL rflu_computegradcellsconstr(pregion,ibegvar,iendvar,ibeggrad, &
1663  iendgrad,varinfo,var,grad)
1664  END IF ! pRegion%grid%nCellsConstr
1665  CASE ( 3 )
1666 ! CALL RFLU_ComputeGradCells_3D(pRegion,iBegVar,iEndVar,iBegGrad, &
1667 ! iEndGrad,var,grad)
1668  CALL rflu_computegradcellsfast_3d(pregion,ibegvar,iendvar,ibeggrad, &
1669  iendgrad,var,grad)
1670 
1671  IF ( pregion%grid%nCellsConstr > 0 ) THEN
1672  CALL rflu_computegradcellsconstr(pregion,ibegvar,iendvar,ibeggrad, &
1673  iendgrad,varinfo,var,grad)
1674  END IF ! pRegion%grid%nCellsConstr
1675  CASE default
1676  CALL errorstop(global,err_reached_default,__line__)
1677  END SELECT ! pMixtInput%stencilDimensCells
1678 
1679 ! ******************************************************************************
1680 ! End
1681 ! ******************************************************************************
1682 
1683  CALL deregisterfunction(global)
1684 
1685  END SUBROUTINE rflu_computegradcellswrapper
1686 
1687 
1688 
1689 
1690 
1691 ! ******************************************************************************
1692 ! End
1693 ! ******************************************************************************
1694 
1696 
1697 
1698 ! ******************************************************************************
1699 !
1700 ! RCS Revision history:
1701 !
1702 ! $Log: RFLU_ModDifferentiationCells.F90,v $
1703 ! Revision 1.12 2008/12/06 08:44:21 mtcampbe
1704 ! Updated license.
1705 !
1706 ! Revision 1.11 2008/11/19 22:17:32 mtcampbe
1707 ! Added Illinois Open Source License/Copyright
1708 !
1709 ! Revision 1.10 2007/02/27 13:03:28 haselbac
1710 ! Enabled 1d computations
1711 !
1712 ! Revision 1.9 2006/04/19 19:42:31 haselbac
1713 ! Added tuned gradient routines
1714 !
1715 ! Revision 1.8 2006/04/07 15:19:19 haselbac
1716 ! Removed tabs
1717 !
1718 ! Revision 1.7 2006/04/07 14:48:12 haselbac
1719 ! Moved WENO routines into separate module, bug fix for 1D grads
1720 !
1721 ! Revision 1.6 2006/03/09 15:05:17 haselbac
1722 ! Bug fix: Missing ampersand
1723 !
1724 ! Revision 1.5 2006/03/09 14:07:10 haselbac
1725 ! Now use generic 1d weight routine
1726 !
1727 ! Revision 1.4 2006/01/06 22:10:00 haselbac
1728 ! Added wrappers and 1D routines
1729 !
1730 ! Revision 1.3 2005/12/25 15:31:33 haselbac
1731 ! Added user-specified constraint weight
1732 !
1733 ! Revision 1.2 2005/10/27 19:03:25 haselbac
1734 ! Separate constrained-gradient routine for performance reasons
1735 !
1736 ! Revision 1.1 2005/10/05 14:33:44 haselbac
1737 ! Initial revision
1738 !
1739 ! ******************************************************************************
1740 
1741 
1742 
1743 
1744 
1745 
1746 
1747 
1748 
1749 
1750 
1751 
1752 
subroutine, public rflu_computewtsx2c_1d(global, m, nMembs, x, z, w)
subroutine rflu_computegradcellsfast_3d(pRegion, iBegVar, iEndVar, iBegGrad, iEndGrad, var, grad)
NT dx
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
subroutine ainv(ajac, ajacin, det, ndim)
Definition: ainv.f90:53
Size order() const
Degree of the element. 1 for linear and 2 for quadratic.
subroutine rflu_invertmatrixsvd(global, nRows, nCols, a, aInv, sCount)
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
subroutine rflu_computegradcells_3d(pRegion, iBegVar, iEndVar, iBegGrad, iEndGrad, var, grad)
subroutine rflu_computegradcells_1d(pRegion, iBegVar, iEndVar, iBegGrad, iEndGrad, var, grad)
double sqrt(double d)
Definition: double.h:73
subroutine rflu_computegradcells_2d(pRegion, iBegVar, iEndVar, iBegGrad, iEndGrad, var, grad)
INTEGER function, public rflu_getconstrtype(pRegion, pPatch, var, ifl)
subroutine rflu_computegradcellsconstr(pRegion, iBegVar, iEndVar, iBegGrad, iEndGrad, varInfo, var, grad)
RT dz() const
Definition: Direction_3.h:133
NT dy
real(rfreal) function, public rflu_getconstrvalue(pRegion, pPatch, var, ifl)
subroutine, public rflu_computegradcellswrapper(pRegion, iBegVar, iEndVar, iBegGrad, iEndGrad, varInfo, var, grad)
subroutine errorstop(global, errorCode, errorLine, addMessage)
Definition: ModError.F90:483
subroutine deregisterfunction(global)
Definition: ModError.F90:469
subroutine rflu_computegradcellsfast_2d(pRegion, iBegVar, iEndVar, iBegGrad, iEndGrad, var, grad)
RT a() const
Definition: Line_2.h:140