Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLU_ModDifferentiationFaces.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 functions at face centroids.
26 !
27 ! Description: None.
28 !
29 ! Notes: None.
30 !
31 ! ******************************************************************************
32 !
33 ! $Id: RFLU_ModDifferentiationFaces.F90,v 1.9 2008/12/06 08:44:21 mtcampbe Exp $
34 !
35 ! Copyright: (c) 2005 by the University of Illinois
36 !
37 ! ******************************************************************************
38 
39 MODULE rflu_moddifferentiationfaces
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 
50  USE rflu_modconstraintutils
51 
52  IMPLICIT NONE
53 
54  PRIVATE
55  PUBLIC :: rflu_computegradfaceswrapper, &
56  rflu_computegradfacesconstr, &
57  rflu_computegradconstrained
58 
59 ! ******************************************************************************
60 ! Declarations and definitions
61 ! ******************************************************************************
62 
63  CHARACTER(CHRLEN) :: RCSIdentString = &
64  '$RCSfile: RFLU_ModDifferentiationFaces.F90,v $ $Revision: 1.9 $'
65 
66 ! ******************************************************************************
67 ! Routines
68 ! ******************************************************************************
69 
70  CONTAINS
71 
72 
73 
74 
75 
76 
77 
78 ! ******************************************************************************
79 !
80 ! Purpose: Compute gradients of any vector or scalar at face centers.
81 !
82 ! Description: None.
83 !
84 ! Input:
85 ! pRegion Pointer to region data
86 ! iBegVar Beginning index of data in var
87 ! iEndVar Ending index of data in var
88 ! iBegGrad Beginning index of data in grad
89 ! iEndGrad Ending index of data in grad
90 ! var Variables of which gradients are to be determined
91 !
92 ! Output:
93 ! grad Gradients of variables at face centers
94 !
95 ! Notes:
96 ! 1. The face gradients differ from the cell gradients in that they are
97 ! computed as weighted sums of variables rather than variable differences
98 ! because there are no variables located at faces.
99 ! 2. If the weighting is changed from inverse-distance to none, then the
100 ! routine RFLU_ComputeStencilMomentsX in RFLU_ModWeights must also be
101 ! adapted.
102 !
103 ! ******************************************************************************
104 
105  SUBROUTINE rflu_computegradfaces(pRegion,iBegVar,iEndVar,iBegGrad,iEndGrad, &
106  var,grad)
107 
108  IMPLICIT NONE
109 
110 ! ******************************************************************************
111 ! Definitions and declarations
112 ! ******************************************************************************
113 
114 ! ==============================================================================
115 ! Arguments
116 ! ==============================================================================
117 
118  INTEGER, INTENT(IN) :: ibegvar,iendvar,ibeggrad,iendgrad
119  REAL(RFREAL), DIMENSION(:,:), POINTER :: var
120  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
121  TYPE(t_region), POINTER :: pregion
122 
123 ! ==============================================================================
124 ! Locals
125 ! ==============================================================================
126 
127  INTEGER :: errorflag,icg,ifg,ifl,igrad,ipatch,isl,ivar
128  REAL(RFREAL) :: c11,c12,c13,c14,c22,c23,c24,c33,c34,c44,dx,dy,dz,r11, &
129  r12,r13,r14,r22,r23,r24,r33,r34,r44,term,term1,term2, &
130  term3,term4,wx,wy,wz
131  REAL(RFREAL) :: fc(xcoord:zcoord)
132  TYPE(t_global), POINTER :: global
133  TYPE(t_grid), POINTER :: pgrid
134  TYPE(t_patch), POINTER :: ppatch
135 
136 ! ******************************************************************************
137 ! Start
138 ! ******************************************************************************
139 
140  global => pregion%global
141 
142  CALL registerfunction(global,'RFLU_ComputeGradFaces',&
143  'RFLU_ModDifferentiationFaces.F90')
144 
145 #ifdef ROCPROF
146  CALL fprofiler_begins("RFLU::ComputeGradFaces")
147 #endif
148 
149  IF ( (iendvar - ibegvar) /= (iendgrad - ibeggrad) ) THEN
150  CALL errorstop(global,err_grad_mismatch,__line__)
151  END IF ! iEndVar
152 
153 ! ******************************************************************************
154 ! Set pointers and variables
155 ! ******************************************************************************
156 
157  pgrid => pregion%grid
158 
159 ! ******************************************************************************
160 ! Loop over faces and compute gradients
161 ! ******************************************************************************
162 
163 ! DEBUG
164 ! DO icg = 1,pGrid%nCellsTot
165 ! DO iVar = iBegVar,iEndVar
166 ! var(iVar,icg) = &
167 ! + REAL(4*(iVar-1) ,RFREAL) &
168 ! + REAL(4*(iVar-1)+1,RFREAL)*pGrid%cofg(XCOORD,icg) &
169 ! + REAL(4*(iVar-1)+2,RFREAL)*pGrid%cofg(YCOORD,icg) &
170 ! + REAL(4*(iVar-1)+3,RFREAL)*pGrid%cofg(ZCOORD,icg)
171 ! END DO ! iVar
172 ! END DO ! icg
173 ! END DEBUG
174 
175 ! ==============================================================================
176 ! Select appropriate dimensionality
177 ! ==============================================================================
178 
179  SELECT CASE ( pregion%mixtInput%dimens )
180 
181 ! ------------------------------------------------------------------------------
182 ! Two dimensions
183 ! ------------------------------------------------------------------------------
184 
185  CASE ( 2 )
186  DO ifg = 1,pgrid%nFaces
187  r11 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_11)
188  r12 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_12)
189  r22 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_22)
190  r13 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_13)
191  r23 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_23)
192  r33 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_33)
193 
194  c11 = 1.0_rfreal/r11
195  c22 = 1.0_rfreal/r22
196  c33 = 1.0_rfreal/r33
197 
198  c12 = - c11*r12
199  c13 = -(c11*r13 + c12*c22*r23)
200 
201  c23 = - c22*r23
202 
203  fc(xcoord) = pgrid%fc(xcoord,ifg)
204  fc(ycoord) = pgrid%fc(ycoord,ifg)
205 
206  DO igrad = ibeggrad,iendgrad
207  grad(xcoord,igrad,ifg) = 0.0_rfreal
208  grad(ycoord,igrad,ifg) = 0.0_rfreal
209  grad(zcoord,igrad,ifg) = 0.0_rfreal
210  END DO ! iVar
211 
212  DO isl = 1,pgrid%f2cs(ifg)%nCellMembs
213  icg = pgrid%f2cs(ifg)%cellMembs(isl)
214 
215  dx = pgrid%cofg(xcoord,icg) - fc(xcoord)
216  dy = pgrid%cofg(ycoord,icg) - fc(ycoord)
217 
218  term = 1.0_rfreal/sqrt(dx*dx + dy*dy)
219 
220  dx = term*dx
221  dy = term*dy
222 
223  term1 = c11*c11*( dx)
224  term2 = c22*c22*( dy + c12*dx)
225  term3 = c33*c33*(term + c23*dy + c13*dx)
226 
227  wx = term*(term1 + c12*term2 + c13*term3)
228  wy = term*( term2 + c23*term3)
229 
230  igrad = ibeggrad
231 
232  DO ivar = ibegvar,iendvar
233  grad(xcoord,igrad,ifg) = grad(xcoord,igrad,ifg) + wx*var(ivar,icg)
234  grad(ycoord,igrad,ifg) = grad(ycoord,igrad,ifg) + wy*var(ivar,icg)
235 
236  igrad = igrad + 1
237  END DO ! iVar
238  END DO ! isl
239  END DO ! ifg
240 
241 ! ------------------------------------------------------------------------------
242 ! Three dimensions
243 ! ------------------------------------------------------------------------------
244 
245  CASE ( 3 )
246  DO ifg = 1,pgrid%nFaces
247  r11 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_11)
248  r12 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_12)
249  r22 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_22)
250  r13 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_13)
251  r23 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_23)
252  r33 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_33)
253  r14 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_14)
254  r24 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_24)
255  r34 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_34)
256  r44 = pgrid%f2cs(ifg)%xyzMoms(xyz_mom_44)
257 
258  c11 = 1.0_rfreal/r11
259  c22 = 1.0_rfreal/r22
260  c33 = 1.0_rfreal/r33
261  c44 = 1.0_rfreal/r44
262 
263  c12 = - c11*r12
264  c13 = -(c11*r13 + c12*c22*r23)
265  c14 = -(c11*r14 + c12*c22*r24 + c13*c33*r34)
266 
267  c23 = - c22*r23
268  c24 = -(c22*r24 + c23*c33*r34)
269 
270  c34 = - c33*r34
271 
272  fc(xcoord) = pgrid%fc(xcoord,ifg)
273  fc(ycoord) = pgrid%fc(ycoord,ifg)
274  fc(zcoord) = pgrid%fc(zcoord,ifg)
275 
276  DO igrad = ibeggrad,iendgrad
277  grad(xcoord,igrad,ifg) = 0.0_rfreal
278  grad(ycoord,igrad,ifg) = 0.0_rfreal
279  grad(zcoord,igrad,ifg) = 0.0_rfreal
280  END DO ! iVar
281 
282  DO isl = 1,pgrid%f2cs(ifg)%nCellMembs
283  icg = pgrid%f2cs(ifg)%cellMembs(isl)
284 
285  dx = pgrid%cofg(xcoord,icg) - fc(xcoord)
286  dy = pgrid%cofg(ycoord,icg) - fc(ycoord)
287  dz = pgrid%cofg(zcoord,icg) - fc(zcoord)
288 
289  term = 1.0_rfreal/sqrt(dx*dx + dy*dy + dz*dz)
290 
291  dx = term*dx
292  dy = term*dy
293  dz = term*dz
294 
295  term1 = c11*c11*( dx)
296  term2 = c22*c22*( dy + c12*dx)
297  term3 = c33*c33*( dz + c23*dy + c13*dx)
298  term4 = c44*c44*(term + c34*dz + c24*dy + c14*dx)
299 
300  wx = term*(term1 + c12*term2 + c13*term3 + c14*term4)
301  wy = term*( term2 + c23*term3 + c24*term4)
302  wz = term*( term3 + c34*term4)
303 
304  igrad = ibeggrad
305 
306  DO ivar = ibegvar,iendvar
307  grad(xcoord,igrad,ifg) = grad(xcoord,igrad,ifg) + wx*var(ivar,icg)
308  grad(ycoord,igrad,ifg) = grad(ycoord,igrad,ifg) + wy*var(ivar,icg)
309  grad(zcoord,igrad,ifg) = grad(zcoord,igrad,ifg) + wz*var(ivar,icg)
310 
311  igrad = igrad + 1
312  END DO ! iVar
313  END DO ! isl
314  END DO ! ifg
315 
316 ! ------------------------------------------------------------------------------
317 ! Default
318 ! ------------------------------------------------------------------------------
319 
320  CASE default
321  CALL errorstop(global,err_reached_default,__line__)
322  END SELECT ! pRegion%mixtInput%dimens
323 
324 ! DEBUG
325 ! DO iGrad = iBegGrad,iEndGrad
326 ! WRITE(*,*) iGrad,MINVAL(grad(XCOORD,iGrad,1:pGrid%nFaces)), &
327 ! MINVAL(grad(YCOORD,iGrad,1:pGrid%nFaces)), &
328 ! MINVAL(grad(ZCOORD,iGrad,1:pGrid%nFaces)), &
329 ! MAXVAL(grad(XCOORD,iGrad,1:pGrid%nFaces)), &
330 ! MAXVAL(grad(YCOORD,iGrad,1:pGrid%nFaces)), &
331 ! MAXVAL(grad(ZCOORD,iGrad,1:pGrid%nFaces))
332 ! END DO ! iGrad
333 ! END DEBUG
334 
335 ! ******************************************************************************
336 ! End
337 ! ******************************************************************************
338 
339 #ifdef ROCPROF
340  CALL fprofiler_ends("RFLU::ComputeGradFaces")
341 #endif
342 
343  CALL deregisterfunction(global)
344 
345  END SUBROUTINE rflu_computegradfaces
346 
347 
348 
349 
350 
351 
352 
353 
354 ! ******************************************************************************
355 !
356 ! Purpose: Compute constrained gradients of any vector or scalar at face
357 ! centers.
358 !
359 ! Description: None.
360 !
361 ! Input:
362 ! pRegion Pointer to region data
363 ! iBegVar Beginning index of data in var
364 ! iEndVar Ending index of data in var
365 ! iBegGrad Beginning index of data in grad
366 ! iEndGrad Ending index of data in grad
367 ! varInfo Variable information
368 ! var Variables of which gradients are to be determined
369 !
370 ! Output:
371 ! grad Gradients of variables at face centers
372 !
373 ! Notes:
374 ! 1. If the weighting is changed from inverse-distance to none, then the
375 ! routine RFLU_ComputeStencilMomentsX in RFLU_ModWeights must also be
376 ! adapted.
377 ! 2. Restricted to linear reconstruction for now.
378 ! 3. Restricted to Dirichlet constraints for now.
379 !
380 ! ******************************************************************************
381 
382  SUBROUTINE rflu_computegradfacesconstr(pRegion,iBegVar,iEndVar,iBegGrad, &
383  iendgrad,varinfo,var,grad)
384 
385  USE modinterfaces, ONLY: rflu_invertmatrixsvd
386 
387  IMPLICIT NONE
388 
389 ! ******************************************************************************
390 ! Definitions and declarations
391 ! ******************************************************************************
392 
393 ! ==============================================================================
394 ! Arguments
395 ! ==============================================================================
396 
397  INTEGER, INTENT(IN) :: ibegvar,iendvar,ibeggrad,iendgrad
398  INTEGER, INTENT(IN) :: varinfo(ibegvar:iendvar)
399  REAL(RFREAL), DIMENSION(:,:), POINTER :: var
400  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
401  TYPE(t_region), POINTER :: pregion
402 
403 ! ==============================================================================
404 ! Locals
405 ! ==============================================================================
406 
407  INTEGER :: errorflag,icg,icol,ifg,ifl,ifl2,igrad,ipatch,isl,irow,ivar, &
408  ncols,nconstr,nrows,scount
409  INTEGER, DIMENSION(:), ALLOCATABLE :: constrtype
410  REAL(RFREAL) :: cwt,dx,dy,dz,term,varc,wtx,wty,wtz,gx,gy,gz
411  REAL(RFREAL) :: colmax(4)
412  REAL(RFREAL) :: fc(xcoord:zcoord)
413  REAL(RFREAL), DIMENSION(:,:), ALLOCATABLE :: a,ainv
414  TYPE(t_global), POINTER :: global
415  TYPE(t_grid), POINTER :: pgrid
416  TYPE(t_patch), POINTER :: ppatch
417 
418 ! ******************************************************************************
419 ! Start
420 ! ******************************************************************************
421 
422  global => pregion%global
423 
424  CALL registerfunction(global,'RFLU_ComputeGradFacesConstr',&
425  'RFLU_ModDifferentiationFaces.F90')
426 
427 #ifdef ROCPROF
428  CALL fprofiler_begins("RFLU::ComputeGradFacesConstr")
429 #endif
430 
431  IF ( (iendvar - ibegvar) /= (iendgrad - ibeggrad) ) THEN
432  CALL errorstop(global,err_grad_mismatch,__line__)
433  END IF ! iEndVar
434 
435 ! ******************************************************************************
436 ! Set pointers and variables
437 ! ******************************************************************************
438 
439  pgrid => pregion%grid
440 
441  cwt = pregion%mixtInput%cReconstFacesWeight
442 
443 ! DEBUG
444 ! DO icg = 1,pGrid%nCellsTot
445 ! DO iVar = iBegVar,iEndVar
446 ! var(iVar,icg) = &
447 ! + REAL(4*(iVar-1) ,RFREAL) &
448 ! + REAL(4*(iVar-1)+1,RFREAL)*pGrid%cofg(XCOORD,icg) &
449 ! + REAL(4*(iVar-1)+2,RFREAL)*pGrid%cofg(YCOORD,icg) &
450 ! + REAL(4*(iVar-1)+3,RFREAL)*pGrid%cofg(ZCOORD,icg)
451 ! END DO ! iVar
452 ! END DO ! icg
453 ! END DEBUG
454 
455 ! ******************************************************************************
456 ! Loop over faces and compute gradients
457 ! ******************************************************************************
458 
459  DO ifl = 1,pgrid%nFacesConstr
460  ifg = pgrid%ifgConstr(ifl)
461 
462 ! ==============================================================================
463 ! Initialize gradients
464 ! ==============================================================================
465 
466  DO igrad = ibeggrad,iendgrad
467  grad(xcoord,igrad,ifg) = 0.0_rfreal
468  grad(ycoord,igrad,ifg) = 0.0_rfreal
469  grad(zcoord,igrad,ifg) = 0.0_rfreal
470  END DO ! iGrad
471 
472  fc(xcoord) = pgrid%fc(xcoord,ifg)
473  fc(ycoord) = pgrid%fc(ycoord,ifg)
474  fc(zcoord) = pgrid%fc(zcoord,ifg)
475 
476  ALLOCATE(constrtype(pgrid%f2cs(ifg)%nBFaceMembs),stat=errorflag)
477  global%error = errorflag
478  IF ( global%error /= err_none ) THEN
479  CALL errorstop(global,err_allocate,__line__,'constrType')
480  END IF ! global%error
481 
482 ! ==============================================================================
483 ! Compute gradients
484 ! ==============================================================================
485 
486  igrad = ibeggrad
487 
488  DO ivar = ibegvar,iendvar
489 
490 ! ------------------------------------------------------------------------------
491 ! Determine number of constraints
492 ! ------------------------------------------------------------------------------
493 
494  nconstr = 0
495 
496  DO isl = 1,pgrid%f2cs(ifg)%nBFaceMembs
497  ipatch = pgrid%f2cs(ifg)%bFaceMembs(1,isl)
498  ifl2 = pgrid%f2cs(ifg)%bFaceMembs(2,isl)
499 
500  ppatch => pregion%patches(ipatch)
501 
502  constrtype(isl) = rflu_getconstrtype(pregion,ppatch,varinfo(ivar), &
503  ifl2)
504 
505  IF ( constrtype(isl) == constr_type_dirichlet ) THEN
506  nconstr = nconstr + 1
507  ELSE
508  constrtype(isl) = constr_type_none
509  END IF ! constrType
510  END DO ! pGrid%f2cs(ifg)%nBFaceMembs
511 
512 ! ------------------------------------------------------------------------------
513 ! Gradients constrained by Dirichlet boundary conditions. Treated as
514 ! soft constraints. NOTE do not need to treat case of unconstrained
515 ! gradients here because they should already have been computed. If
516 ! they have not been computed at this stage, they will be set to zero
517 ! by initialization above.
518 ! ------------------------------------------------------------------------------
519 
520  IF ( nconstr > 0 ) THEN
521 
522 ! ------- Allocate temporary memory --------------------------------------------
523 
524  nrows = pgrid%f2cs(ifg)%nCellMembs + nconstr
525  ncols = pregion%mixtInput%dimens + 1
526 
527  ALLOCATE(a(nrows,ncols),stat=errorflag)
528  global%error = errorflag
529  IF ( global%error /= err_none ) THEN
530  CALL errorstop(global,err_allocate,__line__,'a')
531  END IF ! global%error
532 
533  ALLOCATE(ainv(ncols,nrows),stat=errorflag)
534  global%error = errorflag
535  IF ( global%error /= err_none ) THEN
536  CALL errorstop(global,err_allocate,__line__,'aInv')
537  END IF ! global%error
538 
539 ! ------- Define left-hand side matrix -----------------------------------------
540 
541  SELECT CASE ( pregion%mixtInput%dimens )
542 
543 ! --------- Two dimensions
544 
545  CASE ( 2 )
546  DO isl = 1,pgrid%f2cs(ifg)%nCellMembs
547  icg = pgrid%f2cs(ifg)%cellMembs(isl)
548 
549  dx = pgrid%cofg(xcoord,icg) - fc(xcoord)
550  dy = pgrid%cofg(ycoord,icg) - fc(ycoord)
551 
552  term = 1.0_rfreal/sqrt(dx*dx + dy*dy)
553 
554  a(isl,1) = term
555  a(isl,2) = term*dx
556  a(isl,3) = term*dy
557  END DO ! isl
558 
559  irow = pgrid%f2cs(ifg)%nCellMembs
560 
561  DO isl = 1,pgrid%f2cs(ifg)%nBFaceMembs
562  IF ( constrtype(isl) == constr_type_dirichlet ) THEN
563  ipatch = pgrid%f2cs(ifg)%bFaceMembs(1,isl)
564  ifl2 = pgrid%f2cs(ifg)%bFaceMembs(2,isl)
565 
566  ppatch => pregion%patches(ipatch)
567 
568  dx = ppatch%fc(xcoord,ifl2) - fc(xcoord)
569  dy = ppatch%fc(ycoord,ifl2) - fc(ycoord)
570 
571  term = cwt/sqrt(dx*dx + dy*dy)
572 
573  irow = irow + 1
574 
575  a(irow,1) = term
576  a(irow,2) = term*dx
577  a(irow,3) = term*dy
578  END IF ! constrType
579  END DO ! isl
580 
581 ! --------- Three dimensions
582 
583  CASE ( 3 )
584  DO isl = 1,pgrid%f2cs(ifg)%nCellMembs
585  icg = pgrid%f2cs(ifg)%cellMembs(isl)
586 
587  dx = pgrid%cofg(xcoord,icg) - fc(xcoord)
588  dy = pgrid%cofg(ycoord,icg) - fc(ycoord)
589  dz = pgrid%cofg(zcoord,icg) - fc(zcoord)
590 
591  term = 1.0_rfreal/sqrt(dx*dx + dy*dy + dz*dz)
592 
593  a(isl,1) = term
594  a(isl,2) = term*dx
595  a(isl,3) = term*dy
596  a(isl,4) = term*dz
597  END DO ! isl
598 
599  irow = pgrid%f2cs(ifg)%nCellMembs
600 
601  DO isl = 1,pgrid%f2cs(ifg)%nBFaceMembs
602  IF ( constrtype(isl) == constr_type_dirichlet ) THEN
603  ipatch = pgrid%f2cs(ifg)%bFaceMembs(1,isl)
604  ifl2 = pgrid%f2cs(ifg)%bFaceMembs(2,isl)
605 
606  ppatch => pregion%patches(ipatch)
607 
608  dx = ppatch%fc(xcoord,ifl2) - fc(xcoord)
609  dy = ppatch%fc(ycoord,ifl2) - fc(ycoord)
610  dz = ppatch%fc(zcoord,ifl2) - fc(zcoord)
611 
612  term = cwt/sqrt(dx*dx + dy*dy + dz*dz)
613 
614  irow = irow + 1
615 
616  a(irow,1) = term
617  a(irow,2) = term*dx
618  a(irow,3) = term*dy
619  a(irow,4) = term*dz
620  END IF ! constrType
621  END DO ! isl
622  CASE default
623  CALL errorstop(global,err_reached_default,__line__)
624  END SELECT ! pRegion%mixtInput%dimens
625 
626 ! ------- Compute constrained gradient weights ---------------------------------
627 
628  DO icol = 1,ncols
629  colmax(icol) = -huge(1.0_rfreal)
630 
631  DO irow = 1,nrows
632  colmax(icol) = max(colmax(icol),abs(a(irow,icol)))
633  END DO ! iRow
634 
635  DO irow = 1,nrows
636  a(irow,icol) = a(irow,icol)/colmax(icol)
637  END DO ! iRow
638  END DO ! iCol
639 
640  CALL rflu_invertmatrixsvd(global,nrows,ncols,a,ainv,scount)
641 
642  DO icol = 1,ncols
643  DO irow = 1,nrows
644  ainv(icol,irow) = ainv(icol,irow)/colmax(icol)
645  END DO ! iRow
646  END DO ! iCol
647 
648 ! TEMPORARY
649  IF ( scount /= 0 ) THEN
650  WRITE(*,*) 'ERROR - Singular matrix in RFLU_ComputeGradFacesConstr!'
651  stop
652  END IF ! sCount
653 ! END TEMPORARY
654 
655 ! ------- Compute constrained gradient weights ---------------------------------
656 
657  SELECT CASE ( pregion%mixtInput%dimens )
658  CASE ( 2 )
659  gx = 0.0_rfreal
660  gy = 0.0_rfreal
661 
662  DO isl = 1,pgrid%f2cs(ifg)%nCellMembs
663  icg = pgrid%f2cs(ifg)%cellMembs(isl)
664 
665  dx = pgrid%cofg(xcoord,icg) - fc(xcoord)
666  dy = pgrid%cofg(ycoord,icg) - fc(ycoord)
667 
668  term = 1.0_rfreal/sqrt(dx*dx + dy*dy)
669 
670  gx = gx + term*ainv(2,isl)*var(ivar,icg)
671  gy = gy + term*ainv(3,isl)*var(ivar,icg)
672  END DO ! isl
673 
674  irow = pgrid%f2cs(ifg)%nCellMembs
675 
676  DO isl = 1,pgrid%f2cs(ifg)%nBFaceMembs
677  IF ( constrtype(isl) == constr_type_dirichlet ) THEN
678  ipatch = pgrid%f2cs(ifg)%bFaceMembs(1,isl)
679  ifl2 = pgrid%f2cs(ifg)%bFaceMembs(2,isl)
680 
681  ppatch => pregion%patches(ipatch)
682 
683  dx = ppatch%fc(xcoord,ifl2) - fc(xcoord)
684  dy = ppatch%fc(ycoord,ifl2) - fc(ycoord)
685 
686  term = cwt/sqrt(dx*dx + dy*dy)
687 
688  varc = rflu_getconstrvalue(pregion,ppatch,varinfo(ivar),ifl2)
689 
690  irow = irow + 1
691 
692  gx = gx + term*ainv(2,irow)*varc
693  gy = gy + term*ainv(3,irow)*varc
694  END IF ! constrType
695  END DO ! isl
696 
697  grad(xcoord,igrad,ifg) = gx
698  grad(ycoord,igrad,ifg) = gy
699  grad(zcoord,igrad,ifg) = 0.0_rfreal
700  CASE ( 3 )
701  gx = 0.0_rfreal
702  gy = 0.0_rfreal
703  gz = 0.0_rfreal
704 
705  DO isl = 1,pgrid%f2cs(ifg)%nCellMembs
706  icg = pgrid%f2cs(ifg)%cellMembs(isl)
707 
708  dx = pgrid%cofg(xcoord,icg) - fc(xcoord)
709  dy = pgrid%cofg(ycoord,icg) - fc(ycoord)
710  dz = pgrid%cofg(zcoord,icg) - fc(zcoord)
711 
712  term = 1.0_rfreal/sqrt(dx*dx + dy*dy + dz*dz)
713 
714  gx = gx + term*ainv(2,isl)*var(ivar,icg)
715  gy = gy + term*ainv(3,isl)*var(ivar,icg)
716  gz = gz + term*ainv(4,isl)*var(ivar,icg)
717  END DO ! isl
718 
719  irow = pgrid%f2cs(ifg)%nCellMembs
720 
721  DO isl = 1,pgrid%f2cs(ifg)%nBFaceMembs
722  IF ( constrtype(isl) == constr_type_dirichlet ) THEN
723  ipatch = pgrid%f2cs(ifg)%bFaceMembs(1,isl)
724  ifl2 = pgrid%f2cs(ifg)%bFaceMembs(2,isl)
725 
726  ppatch => pregion%patches(ipatch)
727 
728  dx = ppatch%fc(xcoord,ifl2) - fc(xcoord)
729  dy = ppatch%fc(ycoord,ifl2) - fc(ycoord)
730  dz = ppatch%fc(zcoord,ifl2) - fc(zcoord)
731 
732  term = cwt/sqrt(dx*dx + dy*dy + dz*dz)
733 
734  varc = rflu_getconstrvalue(pregion,ppatch,varinfo(ivar),ifl2)
735 
736  irow = irow + 1
737 
738  gx = gx + term*ainv(2,irow)*varc
739  gy = gy + term*ainv(3,irow)*varc
740  gz = gz + term*ainv(4,irow)*varc
741  END IF ! constrType
742  END DO ! isl
743 
744  grad(xcoord,igrad,ifg) = gx
745  grad(ycoord,igrad,ifg) = gy
746  grad(zcoord,igrad,ifg) = gz
747  CASE default
748  CALL errorstop(global,err_reached_default,__line__)
749  END SELECT ! pRegion%mixtInput%dimens
750 
751 ! ------- Deallocate temporary memory ------------------------------------------
752 
753  DEALLOCATE(a,stat=errorflag)
754  global%error = errorflag
755  IF ( global%error /= err_none ) THEN
756  CALL errorstop(global,err_deallocate,__line__,'a')
757  END IF ! global%error
758 
759  DEALLOCATE(ainv,stat=errorflag)
760  global%error = errorflag
761  IF ( global%error /= err_none ) THEN
762  CALL errorstop(global,err_deallocate,__line__,'aInv')
763  END IF ! global%error
764  END IF ! gradType
765 
766  igrad = igrad + 1
767  END DO ! iVar
768 
769  DEALLOCATE(constrtype,stat=errorflag)
770  global%error = errorflag
771  IF ( global%error /= err_none ) THEN
772  CALL errorstop(global,err_deallocate,__line__,'constrType')
773  END IF ! global%error
774  END DO ! ifl
775 
776 ! ******************************************************************************
777 ! End
778 ! ******************************************************************************
779 
780 ! DEBUG
781 ! DO iGrad = iBegGrad,iEndGrad
782 ! WRITE(*,*) iGrad,MINVAL(grad(XCOORD,iGrad,1:pGrid%nFaces)), &
783 ! MINVAL(grad(YCOORD,iGrad,1:pGrid%nFaces)), &
784 ! MINVAL(grad(ZCOORD,iGrad,1:pGrid%nFaces)), &
785 ! MAXVAL(grad(XCOORD,iGrad,1:pGrid%nFaces)), &
786 ! MAXVAL(grad(YCOORD,iGrad,1:pGrid%nFaces)), &
787 ! MAXVAL(grad(ZCOORD,iGrad,1:pGrid%nFaces))
788 ! END DO ! iGrad
789 ! END DEBUG
790 
791 #ifdef ROCPROF
792  CALL fprofiler_ends("RFLU::ComputeGradFacesConstr")
793 #endif
794 
795  CALL deregisterfunction(global)
796 
797  END SUBROUTINE rflu_computegradfacesconstr
798 
799 
800 
801 
802 
803 
804 
805 ! ******************************************************************************
806 !
807 ! Purpose: Wrapper routine for computing gradients of any vector or scalar at
808 ! face centers.
809 !
810 ! Description: None.
811 !
812 ! Input:
813 ! pRegion Pointer to region data
814 ! iBegVar Beginning index of data in var
815 ! iEndVar Ending index of data in var
816 ! iBegGrad Beginning index of data in grad
817 ! iEndGrad Ending index of data in grad
818 ! var Variables of which gradients are to be determined
819 !
820 ! Output:
821 ! grad Gradients of variables at cell centers
822 !
823 ! Notes: None.
824 !
825 ! ******************************************************************************
826 
827  SUBROUTINE rflu_computegradfaceswrapper(pRegion,iBegVar,iEndVar,iBegGrad, &
828  iendgrad,var,grad)
829 
830  IMPLICIT NONE
831 
832 ! ******************************************************************************
833 ! Definitions and declarations
834 ! ******************************************************************************
835 
836 ! ==============================================================================
837 ! Arguments
838 ! ==============================================================================
839 
840  INTEGER, INTENT(IN) :: ibegvar,iendvar,ibeggrad,iendgrad
841  REAL(RFREAL), DIMENSION(:,:), POINTER :: var
842  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: grad
843  TYPE(t_region), POINTER :: pregion
844 
845 ! ==============================================================================
846 ! Locals
847 ! ==============================================================================
848 
849  TYPE(t_global), POINTER :: global
850  TYPE(t_grid), POINTER :: pgrid
851 
852 ! *****************************************************************************
853 ! Start
854 ! *****************************************************************************
855 
856  global => pregion%global
857 
858  CALL registerfunction(global,'RFLU_ComputeGradFacesWrapper',&
859  'RFLU_ModDifferentiationFaces.F90' )
860 
861 ! ******************************************************************************
862 ! Call gradient routines
863 ! ******************************************************************************
864 
865  SELECT CASE ( pregion%mixtInput%stencilDimensFaces )
866  CASE ( 1 )
867 ! TO DO
868 ! END TO DO
869  CASE ( 2,3 )
870  CALL rflu_computegradfaces(pregion,ibegvar,iendvar,ibeggrad,iendgrad, &
871  var,grad)
872  CASE default
873  CALL errorstop(global,err_reached_default,__line__)
874  END SELECT ! pMixtInput%stencilDimensFaces
875 
876 ! ******************************************************************************
877 ! End
878 ! ******************************************************************************
879 
880  CALL deregisterfunction(global)
881 
882  END SUBROUTINE rflu_computegradfaceswrapper
883 
884 
885 
886 
887 
888 
889 ! ******************************************************************************
890 !
891 ! Purpose: Compute constrained gradients.
892 !
893 ! Description: None.
894 !
895 ! Input:
896 ! global Pointer to global data
897 ! dimens Dimensionality
898 ! nCellMembs Number of cell members in stencil
899 ! nBFaceMembs Number of boundary face members in stencil
900 ! order Order of gradient reconstruction
901 ! dra Coordinate differences of cell members
902 ! drb Coordinate differences of boundary face members
903 ! rhsa Variable values of cell members
904 ! rhsb Variable values of boundary face members
905 !
906 ! Output:
907 ! gradLocal Gradients of variables
908 !
909 ! Notes: None.
910 !
911 ! ******************************************************************************
912 
913  SUBROUTINE rflu_computegradconstrained(global,dimens,nCellMembs,nBFaceMembs, &
914  order,dra,drb,rhsa,rhsb,gradlocal)
915 
916  USE rflu_modstencilsutils, ONLY: rflu_computestencilsize
917  USE modtools, ONLY: compfact
918 
919  IMPLICIT NONE
920 
921 ! ******************************************************************************
922 ! Definitions and declarations
923 ! ******************************************************************************
924 
925 ! ==============================================================================
926 ! Arguments
927 ! ==============================================================================
928 
929  INTEGER, INTENT(IN) :: dimens
930  INTEGER :: nbfacemembs,ncellmembs,order
931  REAL(RFREAL) :: gradlocal(xcoord:xyzmag),rhsa(ncellmembs),rhsb(nbfacemembs)
932  REAL(RFREAL) :: dra(xcoord:zcoord,ncellmembs),drb(xcoord:zcoord,nbfacemembs)
933  TYPE(t_global), POINTER :: global
934 
935 ! ==============================================================================
936 ! Locals
937 ! ==============================================================================
938 
939  INTEGER :: errorflag,isl,j,lda,ldb,ncols,p,q,r,term,workarrayrealsize
940  REAL(RFREAL), DIMENSION(:), ALLOCATABLE :: workarrayreal
941  REAL(RFREAL), DIMENSION(:,:), ALLOCATABLE :: a,b
942 
943 ! ******************************************************************************
944 ! Start
945 ! ******************************************************************************
946 
947  CALL registerfunction(global,'RFLU_ComputeGradConstrained',&
948  'RFLU_ModDifferentiationFaces.F90')
949 
950 ! ******************************************************************************
951 ! Set variables
952 ! ******************************************************************************
953 
954  ncols = rflu_computestencilsize(global,dimens,1,order)
955 
956  workarrayrealsize = nbfacemembs + ncols + 100*ncellmembs
957 
958 ! ******************************************************************************
959 ! Allocate temporary memory
960 ! ******************************************************************************
961 
962  ALLOCATE(a(ncellmembs,ncols),stat=errorflag)
963  global%error = errorflag
964  IF ( global%error /= err_none ) THEN
965  CALL errorstop(global,err_allocate,__line__,'a')
966  END IF ! global%error
967 
968  ALLOCATE(b(nbfacemembs,ncols),stat=errorflag)
969  global%error = errorflag
970  IF ( global%error /= err_none ) THEN
971  CALL errorstop(global,err_allocate,__line__,'b')
972  END IF ! global%error
973 
974  ALLOCATE(workarrayreal(workarrayrealsize),stat=errorflag)
975  global%error = errorflag
976  IF ( global%error /= err_none ) THEN
977  CALL errorstop(global,err_allocate,__line__,'workArrayReal')
978  END IF ! global%error
979 
980 ! ******************************************************************************
981 ! Build matrices
982 ! ******************************************************************************
983 
984 ! ==============================================================================
985 ! Matrix arising from unconstrained part
986 ! ==============================================================================
987 
988  DO isl = 1,ncellmembs
989  a(isl,ncols) = 1.0_rfreal
990 
991  j = 1
992 
993  DO p = 1,order
994  DO q = 0,p
995  DO r = 0,q
996  term = 1.0_rfreal/(compfact(p-q)*compfact(q-r)*compfact(r))
997  a(isl,j) = term*dra(xcoord,isl)**(p-q) &
998  *dra(ycoord,isl)**(q-r) &
999  *dra(zcoord,isl)**r
1000  j = j + 1
1001  END DO ! r
1002  END DO ! q
1003  END DO ! p
1004  END DO ! isl
1005 
1006 ! ==============================================================================
1007 ! Matrix arising from constraints
1008 ! ==============================================================================
1009 
1010  DO isl = 1,nbfacemembs
1011  b(isl,ncols) = 1.0_rfreal
1012 
1013  j = 1
1014 
1015  DO p = 1,order
1016  DO q = 0,p
1017  DO r = 0,q
1018  term = 1.0_rfreal/(compfact(p-q)*compfact(q-r)*compfact(r))
1019  b(isl,j) = term*drb(xcoord,isl)**(p-q) &
1020  *drb(ycoord,isl)**(q-r) &
1021  *drb(zcoord,isl)**r
1022  j = j + 1
1023  END DO ! r
1024  END DO ! q
1025  END DO ! p
1026  END DO ! isl
1027 
1028 ! ******************************************************************************
1029 ! Compute gradients
1030 ! ******************************************************************************
1031 
1032  lda = max(1,ncellmembs)
1033  ldb = max(1,nbfacemembs)
1034 
1035  CALL dgglse(ncellmembs,ncols,nbfacemembs,a,lda,b,ldb,rhsa,rhsb,gradlocal, &
1036  workarrayreal,workarrayrealsize,errorflag)
1037  global%error = errorflag
1038  IF ( global%error /= err_none ) THEN
1039  CALL errorstop(global,err_lapack_output,__line__)
1040  END IF ! global%error
1041 
1042 ! ******************************************************************************
1043 ! Deallocate temporary memory
1044 ! ******************************************************************************
1045 
1046  DEALLOCATE(a,stat=errorflag)
1047  global%error = errorflag
1048  IF ( global%error /= err_none ) THEN
1049  CALL errorstop(global,err_deallocate,__line__,'a')
1050  END IF ! global%error
1051 
1052  DEALLOCATE(b,stat=errorflag)
1053  global%error = errorflag
1054  IF ( global%error /= err_none ) THEN
1055  CALL errorstop(global,err_deallocate,__line__,'b')
1056  END IF ! global%error
1057 
1058  DEALLOCATE(workarrayreal,stat=errorflag)
1059  global%error = errorflag
1060  IF ( global%error /= err_none ) THEN
1061  CALL errorstop(global,err_deallocate,__line__,'workArrayReal')
1062  END IF ! global%error
1063 
1064 ! ******************************************************************************
1065 ! End
1066 ! ******************************************************************************
1067 
1068  CALL deregisterfunction(global)
1069 
1070  END SUBROUTINE rflu_computegradconstrained
1071 
1072 
1073 
1074 
1075 
1076 
1077 ! ******************************************************************************
1078 ! End
1079 ! ******************************************************************************
1080 
1081 END MODULE rflu_moddifferentiationfaces
1082 
1083 
1084 ! ******************************************************************************
1085 !
1086 ! RCS Revision history:
1087 !
1088 ! $Log: RFLU_ModDifferentiationFaces.F90,v $
1089 ! Revision 1.9 2008/12/06 08:44:21 mtcampbe
1090 ! Updated license.
1091 !
1092 ! Revision 1.8 2008/11/19 22:17:32 mtcampbe
1093 ! Added Illinois Open Source License/Copyright
1094 !
1095 ! Revision 1.7 2006/04/07 15:19:19 haselbac
1096 ! Removed tabs
1097 !
1098 ! Revision 1.6 2006/04/07 14:47:03 haselbac
1099 ! Added wrapper func in anticipation of 1D routines
1100 !
1101 ! Revision 1.5 2005/12/25 15:32:10 haselbac
1102 ! Added user-specified constraint weight
1103 !
1104 ! Revision 1.4 2005/10/27 19:12:38 haselbac
1105 ! Bug fixes, separated constr and unconstr grads, bgrads now in own module
1106 !
1107 ! Revision 1.3 2005/10/18 03:00:26 haselbac
1108 ! Bug fix: Incorrect computation of boundary gradients
1109 !
1110 ! Revision 1.2 2005/10/07 20:02:23 haselbac
1111 ! Bug fix: Incorrect loop limits, lead to core dump for parallel cases
1112 !
1113 ! Revision 1.1 2005/10/05 14:33:44 haselbac
1114 ! Initial revision
1115 !
1116 ! ******************************************************************************
1117 
1118 
1119 
1120 
1121 
1122 
1123 
1124 
1125 
1126 
rflu_moddifferentiationfaces::rflu_computegradfaceswrapper
subroutine, public rflu_computegradfaceswrapper(pRegion, iBegVar, iEndVar, iBegGrad, iEndGrad, var, grad)
Definition: RFLU_ModDifferentiationFaces.F90:827
Color::r
unsigned char r() const
Definition: Color.h:68
modtools::compfact
INTEGER function compfact(n)
Definition: ModTools.F90:156
dx
NT dx
Definition: rational_rotation.h:81
Kdtree_d::ExtPoint::term
void term()
Definition: kdtree_d.h:237
max
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
ainv
subroutine ainv(ajac, ajacin, det, ndim)
Definition: ainv.f90:53
Generic_element_2::order
Size order() const
Degree of the element. 1 for linear and 2 for quadratic.
Definition: Generic_element_2.h:117
rflu_invertmatrixsvd
subroutine rflu_invertmatrixsvd(global, nRows, nCols, a, aInv, sCount)
Definition: RFLU_InvertMatrixSVD.F90:50
moderror::registerfunction
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
Color::b
unsigned char b() const
Definition: Color.h:70
modgrid::t_grid
Definition: ModGrid.F90:56
sqrt
double sqrt(double d)
Definition: double.h:73
rflu_moddifferentiationfaces::rflu_computegradfacesconstr
subroutine, public rflu_computegradfacesconstr(pRegion, iBegVar, iEndVar, iBegGrad, iEndGrad, varInfo, var, grad)
Definition: RFLU_ModDifferentiationFaces.F90:382
rflu_modconstraintutils::rflu_getconstrtype
INTEGER function, public rflu_getconstrtype(pRegion, pPatch, var, ifl)
Definition: RFLU_ModConstraintUtils.F90:92
rflu_computestencilsize
INTEGER function rflu_computestencilsize(global, factor, order)
Definition: RFLU_ComputeStencilSize.F90:50
rflu_moddifferentiationfaces::rflu_computegradconstrained
subroutine, public rflu_computegradconstrained(global, dimens, nCellMembs, nBFaceMembs, order, dra, drb, rhsa, rhsb, gradLocal)
Definition: RFLU_ModDifferentiationFaces.F90:913
rflu_modstencilsutils
Definition: RFLU_ModStencilsUtils.F90:39
rflu_modconstraintutils
Definition: RFLU_ModConstraintUtils.F90:39
Direction_3::dz
RT dz() const
Definition: Direction_3.h:133
j
j indices j
Definition: Indexing.h:6
dy
NT dy
Definition: rational_rotation.h:82
q
NT q
Definition: rational_rotation.h:105
rflu_modconstraintutils::rflu_getconstrvalue
real(rfreal) function, public rflu_getconstrvalue(pRegion, pPatch, var, ifl)
Definition: RFLU_ModConstraintUtils.F90:345
rflu_moddifferentiationfaces
Definition: RFLU_ModDifferentiationFaces.F90:39
modglobal::t_global
Definition: ModGlobal.F90:50
moderror::errorstop
subroutine errorstop(global, errorCode, errorLine, addMessage)
Definition: ModError.F90:483
moderror::deregisterfunction
subroutine deregisterfunction(global)
Definition: ModError.F90:469
rflu_moddifferentiationfaces::rflu_computegradfaces
subroutine rflu_computegradfaces(pRegion, iBegVar, iEndVar, iBegGrad, iEndGrad, var, grad)
Definition: RFLU_ModDifferentiationFaces.F90:105
Line_2::a
RT a() const
Definition: Line_2.h:140
modbndpatch::t_patch
Definition: ModBndPatch.F90:91