PlasCom2  1.0
XPACC Multi-physics simluation application
SATUtil.f90
Go to the documentation of this file.
1 MODULE satutil
2 
3  USE grid
4 
5  IMPLICIT NONE
6 
7 CONTAINS
8 
9  SUBROUTINE farfield( &
10  numDim,bufferSizes,numPointsBuffer,patchNormalDir,patchSizes, &
11  numPointsPatch,numPatchPointsOp,patchPointsOp,gridType,gridMetric,&
12  jacobianDeterminant,bcParams,gasParams,rhoBuffer,rhoVBuffer, &
13  rhoEBuffer,viscousFluxBuffer,numscalar,scalarBuffer, &
14  rhoRHS,rhoVRHS,rhoERHS,scalarRHS, &
15  rhoTarget,rhoVTarget,rhoETarget,scalarTarget)
16 
17  IMPLICIT NONE
18 
19  INTEGER(KIND=4) :: numDim,numscalar,patchNormalDir,gridType
20  INTEGER(KIND=8) :: numPointsBuffer,bufferSizes(numdim)
21  INTEGER(KIND=8) :: patchSizes(numdim),numPointsPatch
22  INTEGER(KIND=8) :: numPatchPointsOp
23  INTEGER(KIND=8) :: patchPointsOp(numpatchpointsop)
24  REAL(KIND=8) :: gridMetric(numdim*numdim*numpointsbuffer)
25  REAL(KIND=8) :: jacobianDeterminant(numpointsbuffer)
26  REAL(KIND=8) :: bcParams(3),gasParams(5)
27  REAL(KIND=8) :: rhoBuffer(numpointsbuffer)
28  REAL(KIND=8) :: rhoVBuffer(numdim*numpointsbuffer)
29  REAL(KIND=8) :: rhoEBuffer(numpointsbuffer)
30  REAL(KIND=8) :: viscousFluxBuffer((numdim+2)*numpointsbuffer)
31  REAL(KIND=8) :: scalarBuffer(numscalar*numpointsbuffer)
32  REAL(KIND=8) :: rhoRHS(numpointsbuffer)
33  REAL(KIND=8) :: rhoVRHS(numdim*numpointsbuffer)
34  REAL(KIND=8) :: rhoERHS(numpointsbuffer)
35  REAL(KIND=8) :: scalarRHS(numscalar*numpointsbuffer)
36  REAL(KIND=8) :: rhoTarget(numpointsbuffer)
37  REAL(KIND=8) :: rhoVTarget(numdim*numpointsbuffer)
38  REAL(KIND=8) :: rhoETarget(numpointsbuffer)
39  REAL(KIND=8) :: scalarTarget(numscalar*numpointsbuffer)
40 
41  ! gasParams(:) = Cref Gamma Cp Re
42  ! bcParams(:) = sigma1 sigma2
43 
44  ! ... Local variables
45  INTEGER(4) :: ND, nAuxVars, jj, iDim, iDir
46  INTEGER(8) :: Nc, l0, iPoint, pointIndex, metricOffset
47  INTEGER(4) :: normDir, sgn ! , gas_dv_model
48  REAL(8) :: dsgn, sndspdref2, invtempref, tempref, gamref, spcGasConst
49  REAL(8) :: XI_X, XI_Y, XI_Z, XI_T, bndry_h,sbpBoundaryWeight
50  REAL(8) :: XI_X_TILDE, XI_Y_TILDE, XI_Z_TILDE, gridJacobian
51  REAL(8) :: rhoFac, penaltyFac, scalarPenaltyFac, scalarDiff
52 
53  REAL(8), POINTER :: XIX(:)
54  ! REAL(KIND=RFREAL), Pointer :: XI_TAU(:,:), cv(:,:), gv(:,:)
55  REAL(8), POINTER :: pnorm_vec(:) ! , rhs(:,:)
56  REAL(8), POINTER :: Lambda(:,:), gI1(:), uB(:), Tinv(:,:), Tmat(:,:)
57  REAL(8), POINTER :: Aprime(:,:), gI2(:), matX(:,:), penalty(:) ,correction(:)
58  REAL(8) :: norm_vec(3), gamma, ucon, ubAux,xh
59  REAL(8) :: sigmaI1, sigmaI2, sigmaAux, Sfactor, ibfac_local
60  REAL(8) :: spd_snd, RE, REInv, SAT_sigmaI2_FF, SAT_sigmaI1_FF
61 
62  LOGICAL :: TRIGGERED
63  triggered = .false.
64 
65 
66  ! ... Problem size
67  nd = numdim
68  nc = numpointsbuffer
69  ! N(1:MAX_ND) = 1; Np(1:MAX_ND) = 1;
70  ! Do j = 1, ND
71  ! N(j) = grid_ie(j) - grid_is(j) + 1
72  ! Np(j) = patch_ie(j) - patch_is(j) + 1
73  ! End Do
74  nauxvars = numscalar
75 
76  ! ... Useful quantities about this boundary
77  normdir = abs(patchnormaldir)
78  sgn = normdir / patchnormaldir
79  dsgn = dble(sgn)
80 
81 
82  ! ... Reference quantities
83  ! gasParams = (sndspdref,gamma,Cp,1/Re)
84  ! MJA changed for new nonDimensionalization
85  !sndspdref2 = gasParams(1)*gasParams(1) ! input%sndspdref * input%sndspdref
86  !invtempref = gasParams(3)/sndspdref2 ! input%Cpref / sndspdref2
87  !tempref = 1.0_8 / invtempref
88  !gamref = gasParams(2)
89  !REinv = gasParams(4)
90  gamref = gasparams(2)
91  spcgasconst= gasparams(3)
92  sndspdref2 = gasparams(5)*gasparams(5)
93  tempref = gasparams(4)
94  invtempref = tempref
95  re = gasparams(1)
96  reinv = 0.0_8
97  IF(re .gt. 0) reinv = 1.0_8/re
98 
99 
100  ! ... BC Constants
101  ! bcParams = (sigma1,sigma2,sbpBoundaryStencilWeight)
102  sat_sigmai1_ff = bcparams(1) ! input%SAT_sigmaI1_FF
103  sat_sigmai2_ff = bcparams(2) ! input%SAT_sigmaI2_FF
104  sbpboundaryweight = bcparams(3)
105 
106  ! gas_dv_model = input%gas_dv_model
107 
108 ! WRITE(*,*) 'NORMAL DIR: ',patchNormalDir
109 ! WRITE(*,*) 'SNDSPDREF2 = ',sndspdref2
110 ! WRITE(*,*) 'tempref = ',tempref
111 ! WRITE(*,*) 'gamref = ',gamref
112 ! WRITE(*,*) 'sigma1 = ',SAT_sigmaI1_FF
113 ! WRITE(*,*) 'sigma2 = ',SAT_sigmaI2_FF
114 ! WRITE(*,*) 'REInv = ',REinv
115 ! WRITE(*,*) 'Bweight = ',sbpBoundaryWeight
116 ! WRITE(*,*) 'JAC = ',jacobianDeterminant
117 
118  ! ... memory allocation
119  ALLOCATE(xix(4))
120  ALLOCATE(gi1(nd+2), ub(nd+2), tmat(nd+2,nd+2), tinv(nd+2,nd+2), penalty(nd+2),correction(nd+2))
121  ALLOCATE(gi2(nd+2), aprime(nd+2,nd+2), lambda(nd+2,nd+2), matx(nd+2,nd+2), pnorm_vec(nd))
122 
123 ! WRITE(*,*) 'Farfield ',sgn,normDir
124 
125  xix(:) = 0.0_8
126  IF(gridtype <= unirect) THEN
127  xix(normdir) = gridmetric(normdir)
128  metricoffset = 0
129  gridjacobian = jacobiandeterminant(1)
130  ELSE IF(gridtype == rectilinear) THEN
131  metricoffset = (normdir-1)*numpointsbuffer
132  ELSE
133  metricoffset = (normdir-1)*numdim*numpointsbuffer
134  ENDIF
135 
136 ! WRITE(*,*) 'gridtype,metricoffset = ',gridType,metricOffset
137 
138  ! ... loop over all of the points
139  DO ipoint = 1, numpatchpointsop
140  l0 = patchpointsop(ipoint) + 1 ! assume the points coming in are from C (0-based)
141  ! ... iblank factor (either 1 or zero)
142  ibfac_local = 1.0_8 ! ibfac(l0)
143 
144  ! ... construct the normalized metrics
145  ! XI_X = MT1(l0,t2Map(normDir,1)); XI_Y = 0.0_8; XI_Z = 0.0_8
146  ! if (ND >= 2) XI_Y = MT1(l0,t2Map(normDir,2))
147  ! if (ND == 3) XI_Z = MT1(l0,t2Map(normDir,3))
148  ! XI_T = XI_TAU(l0,normDir)
149  IF(gridtype == curvilinear ) THEN
150  DO idir = 1,numdim
151  xix(idir) = gridmetric(metricoffset + (idir-1)*numpointsbuffer +l0)
152  END DO
153  gridjacobian = jacobiandeterminant(l0)
154  ELSE IF(gridtype == rectilinear) THEN
155  xix(normdir) = gridmetric(metricoffset + l0)
156  gridjacobian = jacobiandeterminant(l0)
157  ENDIF
158 
159  xi_x = xix(1) ! gridMetric(1)
160  xi_y = xix(2)
161  xi_z = xix(3)
162  xi_t = xix(4)
163 
164  ! IF(TRIGGERED.eqv..FALSE.) THEN
165  ! WRITE(*,*) 'XI_X: ',XI_X
166  ! WRITE(*,*) 'XI_Y: ',XI_Y
167  ! WRITE(*,*) 'XI_Z: ',XI_Z
168  ! ENDIF
169 
170  ! ... multiply by the Jacobian
171  xi_x = xi_x * gridjacobian ! JAC(l0)
172  xi_y = xi_y * gridjacobian ! JAC(l0)
173  xi_z = xi_z * gridjacobian ! JAC(l0)
174  xi_t = xi_t * gridjacobian ! JAC(l0)
175 
176  ! ... normalized metrics
177  xh = sqrt((xi_x*xi_x)+(xi_y*xi_y)+(xi_z*xi_z))
178  bndry_h = 1.0_8 / xh
179  xi_x_tilde = xi_x * bndry_h
180  xi_y_tilde = xi_y * bndry_h
181  xi_z_tilde = xi_z * bndry_h
182 
183 
184  ! ... inviscid penalty parameter
185  sigmai1 = -sat_sigmai1_ff * dsgn
186 
187  ! ... pull off the boundary data
188  ! uB(:) = cv(l0,:)
189  ub(1) = rhobuffer(l0)
190  DO idim = 1, numdim
191  ub(1+idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)
192  END DO
193  ub(numdim+2) = rhoebuffer(l0)
194 
195  ! ... target data
196  ! gI1(:) = cvTarget(l0,:)
197  gi1(1) = rhotarget(l0)
198  DO idim = 1, numdim
199  gi1(1+idim) = rhovtarget(l0+(idim-1)*numpointsbuffer)
200  END DO
201  gi1(numdim+2) = rhoetarget(l0)
202 
203  ! ... follow Magnus's algorithm (norm_vec has unit length)
204  norm_vec(1:3) = (/ xi_x_tilde, xi_y_tilde, xi_z_tilde /)
205  pnorm_vec(:) = norm_vec(1:nd) * xh
206 
207  ! ... use target data to make implicit work better?
208  gamma = gamref ! gv(l0,1)
209  ! If (gas_dv_model == DV_MODEL_IDEALGAS .or. &
210  ! gas_dv_model == DV_MODEL_IDEALGAS_MIXTURE) then
211  ! IF(TRIGGERED.eqv..FALSE.) THEN
212  ! WRITE(*,*) 'XI_X_TILDE: ',XI_X_TILDE
213  ! WRITE(*,*) 'XI_Y_TILDE: ',XI_Y_TILDE
214  ! WRITE(*,*) 'XI_Z_TILDE: ',XI_Z_TILDE
215  ! WRITE(*,*) 'XI_X: ',XI_X
216  ! WRITE(*,*) 'XI_Y: ',XI_Y
217  ! WRITE(*,*) 'XI_Z: ',XI_Z
218  ! WRITE(*,*) 'BOUNDARYH: ',bndry_h
219  ! WRITE(*,*) 'GridJacobian: ',gridJacobian
220  ! WRITE(*,*) 'STATE: ',uB(:)
221  ! WRITE(*,*) 'TARGET:',gI1(:)
222 
223  ! ENDIF
224 
225  CALL sat_form_roe_matrices(nd, ub, gi1, gamma, pnorm_vec, tmat, tinv, lambda)
226 
227  ! Else If (gas_dv_model == DV_MODEL_THERMALLY_PERFECT) then
228  !
229  ! Call SAT_Form_Roe_Matrices_TP(uB, gI1, ND, gamma, gamref, pnorm_vec, Tmat, Tinv, Lambda)
230  ! End If
231 
232  ! ... save ucon, speed_sound
233  ucon = lambda(1,1)
234  spd_snd = lambda(nd+1,nd+1) - ucon
235 
236  ! ... only pick those Lambda that are incoming
237  If ( sgn == 1 ) Then
238  Do jj = 1, nd+2
239  lambda(jj,jj) = max(lambda(jj,jj),0.0_8)
240  End Do
241  Else
242  Do jj = 1, nd+2
243  lambda(jj,jj) = min(lambda(jj,jj),0.0_8)
244  End Do
245  End If
246 
247  ! ... modify Lambda if subsonic outflow
248  ! ... left boundary
249  if (sgn == 1 .and. &
250  ucon < 0.0_8 .and. ucon + spd_snd > 0.0_8) then
251  lambda(nd+2,nd+1) = lambda(nd+2,nd+1) - ucon - spd_snd
252  ! right boundary
253  else if (sgn == -1 .and. &
254  ucon > 0.0_8 .and. ucon - spd_snd < 0.0_8) then
255  lambda(nd+1,nd+2) = lambda(nd+1,nd+2) - ucon + spd_snd
256  end if
257 
258  ! ... multiply Lambda' into X
259  matx = matmul(lambda,tinv)
260 
261  ! ... compute the characteristic matrix
262  aprime = matmul(tmat,matx)
263 
264  ! ... subtract off the target
265  ub(1:nd+2) = ub(1:nd+2) - gi1(1:nd+2)
266 
267  ! ... compute the characteristic penalty vector
268  gi2(1:nd+2) = matmul(aprime,ub)
269 
270  penaltyfac = sbpboundaryweight*sigmai1
271  correction(:) = penaltyfac*gi2(:)
272  ! ... compute penalty (Bndry_h is included in Lambda already)
273  ! Do jj = 1, ND+2
274  ! rhs(l0,jj) = rhs(l0,jj) + penaltyFac * gI2(jj)
275  ! End Do
276  rhorhs(l0) = rhorhs(l0) + correction(1)
277  DO idim = 1,numdim
278  pointindex = l0 + (idim-1)*numpointsbuffer
279  rhovrhs(pointindex) = rhovrhs(pointindex) + correction(1+idim)
280  END DO
281  rhoerhs(l0) = rhoerhs(l0) + correction(numdim+2)
282 
283 ! IF(TRIGGERED.eqv..FALSE.) THEN
284 ! WRITE(*,*) 'RHS: ',rhoRHS(l0),rhoVRHS(l0),rhoVRHS(l0+numPointsBuffer),rhoERHS(l0)
285 ! WRITE(*,*) 'CORRECTION: ',correction(:)
286 ! TRIGGERED = .TRUE.
287 ! ENDIF
288 
289  ! ... SAT FARFIELD (species treatment)
290  IF(numscalar > 0) THEN
291  If (dsgn * ucon > 0) Then
292  sfactor = ucon * ibfac_local * penaltyfac
293  rhofac = rhobuffer(l0)/rhotarget(l0)
294  scalarpenaltyfac = sfactor*rhofac
295  Do jj = 1, numscalar
296  pointindex = (jj-1)*numpointsbuffer + l0
297  scalardiff = (scalarbuffer(pointindex) - scalartarget(pointindex))
298  scalarrhs(pointindex) = scalarrhs(pointindex) + scalarpenaltyfac * scalardiff
299  ! rhs_auxVars(l0,jj) = rhs_auxVars(l0,jj) + Sfactor * &
300  ! (auxVars(l0,jj) - auxVarsTarget(l0,jj)/cvTarget(l0,1)*cv(l0,1))
301  End do
302  End If
303  END IF
304 
305  ! VISCOUS PART
306  IF(reinv > 0.0) THEN
307  ! ... factor
308 ! sigmaI2 = dsgn * SAT_sigmaI2_FF * SBP_invPdiag_block1(1) / bndry_h !* JAC(l0)
309  sigmai2 = dsgn * sat_sigmai2_ff*sbpboundaryweight/bndry_h !* JAC(l0)
310 
311  ! ... target state
312  do jj = 1, nd+2
313  gi2(jj) = 0.0_8
314  end do
315 
316  ! ... boundary data (already includes 1/Re factor)
317  ub(:) = viscousfluxbuffer(l0) * xi_x_tilde
318  if (nd >= 2) &
319  ub(:) = ub(:) + viscousfluxbuffer(l0) * xi_y_tilde
320  if (nd == 3) &
321  ub(:) = ub(:) + viscousfluxbuffer(l0) * xi_z_tilde
322  ub(:) = 0
323  ! ... penalty term
324  penalty(:) = sigmai2 * (ub(:) - gi2(:))
325 
326  ! ... add to the rhs
327  rhorhs(l0) = rhorhs(l0) + ibfac_local * penalty(1)
328  DO idim = 1,numdim
329  pointindex = l0 + (idim-1)*numpointsbuffer
330  rhovrhs(pointindex) = rhovrhs(pointindex) + ibfac_local*penalty(1+idim)
331  END DO
332  rhoerhs(l0) = rhoerhs(l0) + ibfac_local * penalty(numdim+2)
333 
334  ! species transport
335 ! sigmaAux = ibfac_local * sigmaI2
336 ! do jj = 1, nAuxVars
337 ! ubAux = patch_ViscousFluxAux(jj,lp,normDir)
338 ! rhs_auxVars(l0,jj) = rhs_auxVars(l0,jj) + sigmaAux* ubAux
339 ! end do
340 
341  END IF
342  END DO
343 
344  DEALLOCATE(xix)
345  DEALLOCATE(correction)
346  DEALLOCATE(gi1, ub, tmat, tinv, penalty)
347  DEALLOCATE(gi2, aprime, lambda, matx, pnorm_vec)
348 
349  RETURN
350 
351  END SUBROUTINE farfield
352 
353  SUBROUTINE slip_adiabatic( &
354  numDim,bufferSizes,numPointsBuffer,patchNormalDir,patchSizes, &
355  numPointsPatch,numPatchPointsOp,patchPointsOp,gridType, &
356  gridMetric,jacobianDeterminant,bcParams,gasParams, &
357  rhoBuffer,rhoVBuffer,rhoEBuffer,numscalar,scalarBuffer, &
358  rhoRHS,rhoVRHS,rhoERHS, scalarRHS,rhoTarget,rhoVTarget, &
359  rhoETarget,scalarTarget)
360 
361  IMPLICIT NONE
362 
363  INTEGER(KIND=4) :: numDim,numscalar,patchNormalDir,gridType
364  INTEGER(KIND=8) :: numPointsBuffer,bufferSizes(numdim)
365  INTEGER(KIND=8) :: patchSizes(numdim),numPointsPatch
366  INTEGER(KIND=8) :: numPatchPointsOp
367  INTEGER(KIND=8) :: patchPointsOp(numpatchpointsop)
368  REAL(KIND=8) :: gridMetric(numdim*numdim*numpointsbuffer)
369  REAL(KIND=8) :: jacobianDeterminant(numpointsbuffer)
370  REAL(KIND=8) :: bcParams(3),gasParams(5)
371  REAL(KIND=8) :: rhoBuffer(numpointsbuffer)
372  REAL(KIND=8) :: rhoVBuffer(numdim*numpointsbuffer)
373  REAL(KIND=8) :: rhoEBuffer(numpointsbuffer)
374  REAL(KIND=8) :: scalarBuffer(numscalar*numpointsbuffer)
375  REAL(KIND=8) :: rhoRHS(numpointsbuffer)
376  REAL(KIND=8) :: rhoVRHS(numdim*numpointsbuffer)
377  REAL(KIND=8) :: rhoERHS(numpointsbuffer)
378  REAL(KIND=8) :: scalarRHS(numscalar*numpointsbuffer)
379  REAL(KIND=8) :: rhoTarget(numpointsbuffer)
380  REAL(KIND=8) :: rhoVTarget(numdim*numpointsbuffer)
381  REAL(KIND=8) :: rhoETarget(numpointsbuffer)
382  REAL(KIND=8) :: scalarTarget(numscalar*numpointsbuffer)
383 
384 
385  ! ... Local variables
386  INTEGER(4) :: ND, nAuxVars, jj, iDim, iDir
387  INTEGER(8) :: Nc, l0, iPoint, pointIndex,metricOffset
388  INTEGER(4) :: normDir, sgn ! , gas_dv_model
389  ! INTEGER, Dimension(MAX_ND) :: N, Np
390  ! INTEGER, Pointer :: iblank(:), t2Map(:,:)
391  ! INTEGER, Pointer :: grid_is(:), grid_ie(:), patch_is(:), patch_ie(:)
392  ! TYPE (t_grid), Pointer :: grid
393  ! TYPE (t_mixt), Pointer :: state
394  ! TYPE (t_mixt_input), Pointer :: input
395  REAL(8) :: dsgn, sndspdref2, invtempref, tempref, gamref, spcGasConst
396  REAL(8) :: XI_X, XI_Y, XI_Z, XI_T, bndry_h,sbpBoundaryWeight
397  REAL(8) :: XI_X_TILDE, XI_Y_TILDE, XI_Z_TILDE, gridJacobian
398  REAL(8) :: rhoFac, penaltyFac, scalarPenaltyFac, scalarDiff
399 
400  REAL(8), POINTER :: XIX(:), bVelocity(:), vWall(:), vWallTarget(:), vTarget(:)
401  ! REAL(KIND=RFREAL), Pointer :: XI_TAU(:,:), cv(:,:), gv(:,:)
402  REAL(8), POINTER :: pnorm_vec(:) ! , rhs(:,:)
403  ! REAL(KIND=RFREAL), Pointer :: patch_ViscousFlux(:,:,:)
404  ! REAL(KIND=RFREAL), Pointer :: patch_ViscousFluxAux(:,:,:)
405  ! REAL(KIND=RFREAL), Pointer :: patch_ViscousRHSAux(:,:,:), rhs_auxVars(:,:)
406  ! REAL(KIND=RFREAL), Pointer :: auxVars(:,:), Lambda(:,:), ibfac(:)
407  REAL(8), POINTER :: Lambda(:,:), gI1(:), uB(:), Tinv(:,:), Tmat(:,:)
408  ! REAL(KIND=RFREAL), Pointer :: gI1(:), ub(:), Tinv(:,:), Tmat(:,:)
409  REAL(8), POINTER :: Aprime(:,:), gI2(:), matX(:,:), penalty(:) ,correction(:)
410  ! REAL(KIND=RFREAL), Pointer :: auxVarsTarget(:,:), cvTarget(:,:)
411  ! REAL(KIND=RFREAL), Pointer :: SBP_invPdiag_block1(:), JAC(:), MT1(:,:)
412  REAL(8) :: norm_vec(3), gamma, uCon, ubAux, specRad, vDotXi, xh
413  REAL(8) :: sigmaI1, sigmaI2, sigmaAux, Sfactor, ibfac_local
414  REAL(8) :: spd_snd, RE, REInv, SAT_sigmaI1,pressureTarget
415 
416  LOGICAL :: TRIGGERED
417  triggered = .false.
418 
419  ! ... Main data containers
420  ! patch_gridID = patch%gridID
421  ! grid => region%grid(patch_gridID)
422  ! state => region%state(patch_gridID)
423  ! input => grid%input
424 
425  ! ... Pointer dereferences
426  ! grid_is => grid%is
427  ! grid_ie => grid%ie
428  ! patch_is => patch%is
429  ! patch_ie => patch%ie
430  ! cv => state%cv
431  ! gv => state%gv
432  ! cvTarget => state%cvTarget
433  ! rhs => state%rhs
434  ! iblank => grid%iblank
435  ! JAC => grid%JAC
436  ! MT1 => grid%MT1
437  ! ibfac => grid%ibfac
438  ! t2Map => region%global%t2Map
439  ! XI_TAU => grid%XI_TAU
440  ! rhs_auxVars => state%rhs_auxVars
441  ! patch_ViscousFlux => patch%ViscousFlux
442  ! patch_ViscousFluxAux => patch%ViscousFluxAux
443  ! SBP_invPdiag_block1 => grid%SBP_invPdiag_block1
444  ! auxVars => state%auxVars
445  ! auxVarsTarget => state%auxVarsTarget
446 
447  ! ... Problem size
448  nd = numdim
449  nc = numpointsbuffer
450 
451  nauxvars = numscalar
452 
453  ! ... Useful quantities about this boundary
454  normdir = abs(patchnormaldir)
455  sgn = normdir / patchnormaldir
456  dsgn = dble(sgn)
457 
458 
459  ! ... Reference quantities
460  ! gasParams = (sndspdref,gamma,Cp,1/Re)
461  ! bcParams(:) = sigma1 sigma2
462  ! MJA changed for new nonDimensionalization
463  !sndspdref2 = gasParams(1)*gasParams(1) ! input%sndspdref * input%sndspdref
464  !invtempref = gasParams(3)/sndspdref2 ! input%Cpref / sndspdref2
465  !tempref = 1.0_8 / invtempref
466  !gamref = gasParams(2)
467  !REinv = gasParams(4)
468  gamref = gasparams(2)
469  spcgasconst= gasparams(3)
470  sndspdref2 = gasparams(5)*gasparams(5)
471  tempref = gasparams(4)
472  invtempref = tempref
473  re = gasparams(1)
474  reinv = 1.0_8/re
475 
476  ! ... BC Constants
477  ! bcParams = (sigma1,sbpBoundaryStencilWeight)
478  sat_sigmai1 = bcparams(1) ! input%SAT_sigmaI1_FF
479  sbpboundaryweight = bcparams(2)
480 
481  ! gas_dv_model = input%gas_dv_model
482 
483  ! ... memory allocation
484  ALLOCATE(xix(4),vwall(nd),bvelocity(nd),vwalltarget(nd),vtarget(nd))
485  ALLOCATE(gi1(nd+2), ub(nd+2), tmat(nd+2,nd+2), tinv(nd+2,nd+2), penalty(nd+2),correction(nd+2))
486  ALLOCATE(gi2(nd+2), aprime(nd+2,nd+2), lambda(nd+2,nd+2), matx(nd+2,nd+2), pnorm_vec(nd))
487 
488  xix(:) = 0.0_8
489  IF(gridtype <= unirect) THEN
490  xix(normdir) = gridmetric(normdir)
491  metricoffset = 0
492  gridjacobian = jacobiandeterminant(1)
493  ELSE IF(gridtype == rectilinear) THEN
494  metricoffset = (normdir-1)*numpointsbuffer
495  ELSE
496  metricoffset = (normdir-1)*numdim*numpointsbuffer
497  ENDIF
498 
499 ! XIX(:) = 0.0_8
500 ! XIX(normDir) = gridMetric(normDir) ! only diagonal terms are non-zero
501 
502 ! WRITE(*,*) 'XIX = ',XIX
503 ! WRITE(*,*) 'NORMAL DIR: ',patchNormalDir
504 ! WRITE(*,*) 'SNDSPDREF2 = ',sndspdref2
505 ! WRITE(*,*) 'tempref = ',tempref
506 ! WRITE(*,*) 'gamref = ',gamref
507 ! WRITE(*,*) 'sigma1 = ',SAT_sigmaI1
508 ! WRITE(*,*) 'REInv = ',REinv
509 ! WRITE(*,*) 'Bweight = ',sbpBoundaryWeight
510 ! WRITE(*,*) 'JAC = ',jacobianDeterminant
511 
512 ! TRIGGERED = .TRUE.
513 
514  ! ... loop over all of the points
515  DO ipoint = 1, numpatchpointsop
516 
517  l0 = patchpointsop(ipoint) + 1 ! assume the points coming in are from C (0-based)
518  ! lp = (k-patch_is(3))*Np(1)*Np(2) + (j-patch_is(2))*Np(1) + (i-patch_is(1)+1)
519 
520  ! IF(iPoint == (numPatchPointsOp/2)) TRIGGERED = .FALSE.
521 
522  ! ... iblank factor (either 1 or zero)
523  ibfac_local = 1.0_8 ! ibfac(l0)
524 
525  ! ... construct the normalized metrics
526  ! XI_X = MT1(l0,t2Map(normDir,1)); XI_Y = 0.0_8; XI_Z = 0.0_8
527  ! if (ND >= 2) XI_Y = MT1(l0,t2Map(normDir,2))
528  ! if (ND == 3) XI_Z = MT1(l0,t2Map(normDir,3))
529  ! XI_T = XI_TAU(l0,normDir)
530  ! IF(ND > 1) XI_Y = gridMetric(2)
531  ! IF(ND > 2) XI_Z = gridMetric(3)
532 
533 
534  IF(gridtype == curvilinear ) THEN
535  DO idir = 1,numdim
536  xix(idir) = gridmetric(metricoffset + (idir-1)*numpointsbuffer +l0)
537  END DO
538  gridjacobian = jacobiandeterminant(l0)
539  ELSE IF(gridtype == rectilinear) THEN
540  xix(normdir) = gridmetric(metricoffset + l0)
541  gridjacobian = jacobiandeterminant(l0)
542  ENDIF
543 
544  ! ... multiply by the Jacobian
545  xi_x = xix(1) * gridjacobian ! JAC(l0)
546  xi_y = xix(2) * gridjacobian ! JAC(l0)
547  xi_z = xix(3) * gridjacobian ! JAC(l0)
548  xi_t = xix(4) * gridjacobian ! JAC(l0)
549 
550  ! ... normalized metrics
551  xh = sqrt( (xi_x * xi_x) + (xi_y * xi_y) + (xi_z * xi_z) )
552  bndry_h = 1.0_8 / xh
553  xi_x_tilde = xi_x * bndry_h
554  xi_y_tilde = xi_y * bndry_h
555  xi_z_tilde = xi_z * bndry_h
556 
557 
558  ! ... follow Magnus's algorithm (norm_vec has unit length)
559  norm_vec(1:3) = (/ xi_x_tilde, xi_y_tilde, xi_z_tilde /)
560  pnorm_vec(:) = norm_vec(1:nd) * xh
561 
562 
563  ! ... penalty parameter
564  sigmai1 = -sat_sigmai1 * dsgn
565 
566  ! ... inviscid penalty parameter
567  ! sigmaI1 = -SAT_sigmaI1_FF * dsgn
568 
569  ! ... pull off the boundary data
570  ! uB(:) = cv(l0,:)
571  ub(1) = rhobuffer(l0)
572  vdotxi = 0
573  DO idim = 1, numdim
574  ub(1+idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)
575  bvelocity(idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)/rhobuffer(l0)
576  vdotxi = vdotxi + bvelocity(idim)*xix(idim)*gridjacobian*bndry_h
577  END DO
578  ub(numdim+2) = rhoebuffer(l0)
579 
580  vwall(1:numdim) = vdotxi*norm_vec(1:numdim)
581  ! NON MOVING GRIDS!!
582  vwalltarget(1:numdim) = 0.0_8
583  ! = dot_product(XYZ_TAU(l0,1:ND),norm_vec(1:ND)) * norm_vec(1:ND)
584 
585 
586  DO idim = 1,numdim
587  vtarget(idim) = bvelocity(idim) - ( vwall(idim) - vwalltarget(idim) )
588  END DO
589 ! pressureTarget = pressureBuffer(l0)
590 
591  ! ... target data
592  ! gI1(:) = cvTarget(l0,:)
593  gi1(1) = ub(1)
594  gi1(numdim+2) = ub(numdim+2)
595  DO idim = 1, numdim
596  gi1(1+idim) = gi1(1) * vtarget(idim)
597  gi1(numdim+2) = gi1(numdim+2) + 0.5_8*gi1(1)*(vtarget(idim)**2 - bvelocity(idim)**2)
598  END DO
599 
600 
601  ! ... use target data to make implicit work better?
602  gamma = gamref ! gv(l0,1)
603  ! If (gas_dv_model == DV_MODEL_IDEALGAS .or. &
604  ! gas_dv_model == DV_MODEL_IDEALGAS_MIXTURE) then
605 
606 
607 ! IF(TRIGGERED .EQV. .FALSE.) THEN
608 ! WRITE(*,*) 'L0 = ',l0
609 ! WRITE(*,*) 'XIX = ',XIX
610 ! WRITE(*,*) '1/xh = ',1.0_8/xh
611 ! WRITE(*,*) 'VdotXi = ',vdotxi
612 ! WRITE(*,*) 'bVelocity = ',bVelocity
613 ! WRITE(*,*) 'pnorm: ',pnorm_vec
614 ! WRITE(*,*) 'State: ',uB
615 ! WRITE(*,*) 'Target: ',gI1
616 ! TRIGGERED = .TRUE.
617 ! ENDIF
618 
619  CALL sat_form_roe_matrices(nd, ub, gi1, gamma, pnorm_vec, tmat, tinv, lambda)
620 
621  specrad = 0.0_8
622  ucon = vdotxi*xh
623 
624  ! ... only pick those Lambda that are incoming
625  If ( sgn == 1 ) Then
626  Do jj = 1, nd+2
627  lambda(jj,jj) = max(lambda(jj,jj),0.0_8)
628  specrad = max(abs(lambda(jj,jj)),abs(specrad))
629  End Do
630  Else
631  Do jj = 1, nd+2
632  lambda(jj,jj) = min(lambda(jj,jj),0.0_8)
633  specrad = max(abs(lambda(jj,jj)),abs(specrad))
634  End Do
635  End If
636 
637 ! IF(TRIGGERED.eqv..FALSE.) THEN
638 ! WRITE(*,*) 'XI_X_TILDE: ',XI_X_TILDE
639 ! WRITE(*,*) 'XI_Y_TILDE: ',XI_Y_TILDE
640 ! WRITE(*,*) 'XI_Z_TILDE: ',XI_Z_TILDE
641 ! WRITE(*,*) 'XI_X: ',XI_X
642 ! WRITE(*,*) 'XI_Y: ',XI_Y
643 ! WRITE(*,*) 'XI_Z: ',XI_Z
644 ! WRITE(*,*) 'BOUNDARYH: ',bndry_h
645 ! WRITE(*,*) 'STATE: ',uB(:)
646 ! WRITE(*,*) 'TARGET:',gI1(:)
647 ! WRITE(*,*) 'UCON: ',uCon
648 ! WRITE(*,*) 'SpecRad: ',specRad
649 ! WRITE(*,*) 'vDotXi: ',vDotXi
650 ! ENDIF
651 
652  matx = matmul(lambda,tinv)
653 
654  ! ... compute the characteristic matrix
655  aprime = matmul(tmat,matx)
656 
657  ! ... subtract off the target
658  ub(1:nd+2) = ub(1:nd+2) - gi1(1:nd+2)
659 
660  ! ... compute the characteristic penalty vector
661  gi2(1:nd+2) = matmul(aprime,ub)
662 
663  penaltyfac = sbpboundaryweight*sigmai1
664  correction(:) = penaltyfac*gi2(:)
665 
666  ! IBFAC would be considered here
667  rhorhs(l0) = rhorhs(l0) + ibfac_local * correction(1)
668  DO idim = 1,numdim
669  pointindex = l0 + (idim-1)*numpointsbuffer
670  rhovrhs(pointindex) = rhovrhs(pointindex) + ibfac_local*correction(1+idim)
671  END DO
672  rhoerhs(l0) = rhoerhs(l0) + ibfac_local * correction(numdim+2)
673 
674 ! IF(TRIGGERED.eqv..FALSE.) THEN
675 ! WRITE(*,*) 'CORRECTION(',l0,'): ',correction(:)
676 ! TRIGGERED = .TRUE.
677 ! ENDIF
678 
679  ! ... SAT_SLIP_ADIABATIC (species treatment)
680  IF(numscalar > 0) THEN
681  If (dsgn * ucon > 0) Then
682  sfactor = -ucon * dsgn * ibfac_local * sbpboundaryweight
683  Do jj = 1, numscalar
684  pointindex = (jj-1)*numpointsbuffer + l0
685  scalarrhs(pointindex) = scalarrhs(pointindex) + sfactor*scalarbuffer(pointindex)
686  End do
687  End If
688  END IF
689 
691  ! VISCOUS PART OMITTED (for now)
692 
693 
694 ! if ( REinv > 0.0_rfreal ) then
695 ! sigmaI2 = ibfac_local * dsgn * 1_8 * SBP_invPdiag_block1(1) * JAC(l0)
696 ! gI2 = 0_8
697 ! uB(:) = patch_ViscousFlux(:,lp,normDir)
698 ! penalty(:) = sigmaI2 * (uB(:) - gI2(:))
699 ! penalty(normDir+1) = 0_8 !normal momentum
700 ! do jj = 1, ND+2
701 ! rhs(l0,jj) = rhs(l0,jj) + penalty(jj)
702 ! end do
703 ! do jj = 1, nAuxVars
704 ! ubAux = patch_ViscousFluxAux(jj,lp,normDir)
705 ! rhs_auxVars(l0,jj) = rhs_auxVars(l0,jj) + sigmaI2* ubAux
706 ! end do
707 ! endif !REinv > 0.0_rfreal
708 
709 
710  END DO
711 
712  DEALLOCATE(xix,vwall,bvelocity,vwalltarget,vtarget)
713  DEALLOCATE(gi1, ub, tmat, tinv, penalty, correction)
714  DEALLOCATE(gi2, aprime, lambda, matx, pnorm_vec)
715 
716  RETURN
717 
718  END SUBROUTINE slip_adiabatic
719 
720  SUBROUTINE noslip_isothermal( &
721  numDim,bufferSizes,numPointsBuffer,patchNormalDir,patchSizes, &
722  numPointsPatch,numPatchPointsOp,patchPointsOp,gridType, &
723  gridMetric,jacobianDeterminant,bcParams,gasParams, &
724  rhoBuffer,rhoVBuffer,rhoEBuffer,numscalar,scalarBuffer, &
725  rhoRHS,rhoVRHS,rhoERHS, scalarRHS,rhoTarget,rhoVTarget, &
726  rhoETarget,scalarTarget,muBuffer,lambdaBuffer)
727 
728  IMPLICIT NONE
729 
730  INTEGER(KIND=4) :: numDim,numscalar,patchNormalDir,gridType
731  INTEGER(KIND=8) :: numPointsBuffer,bufferSizes(numdim)
732  INTEGER(KIND=8) :: patchSizes(numdim),numPointsPatch
733  INTEGER(KIND=8) :: numPatchPointsOp
734  INTEGER(KIND=8) :: patchPointsOp(numpatchpointsop)
735  REAL(KIND=8) :: gridMetric(numdim*numdim*numpointsbuffer)
736  REAL(KIND=8) :: jacobianDeterminant(numpointsbuffer)
737  REAL(KIND=8) :: bcParams(4),gasParams(6)
738  REAL(KIND=8) :: rhoBuffer(numpointsbuffer)
739  REAL(KIND=8) :: rhoVBuffer(numdim*numpointsbuffer)
740  REAL(KIND=8) :: rhoEBuffer(numpointsbuffer)
741  REAL(KIND=8) :: scalarBuffer(numscalar*numpointsbuffer)
742  REAL(KIND=8) :: rhoRHS(numpointsbuffer)
743  REAL(KIND=8) :: rhoVRHS(numdim*numpointsbuffer)
744  REAL(KIND=8) :: rhoERHS(numpointsbuffer)
745  REAL(KIND=8) :: scalarRHS(numscalar*numpointsbuffer)
746  REAL(KIND=8) :: rhoTarget(numpointsbuffer)
747  REAL(KIND=8) :: rhoVTarget(numdim*numpointsbuffer)
748  REAL(KIND=8) :: rhoETarget(numpointsbuffer)
749  REAL(KIND=8) :: scalarTarget(numscalar*numpointsbuffer)
750  REAL(KIND=8) :: muBuffer(numpointsbuffer)
751  REAL(KIND=8) :: lambdaBuffer(numpointsbuffer)
752 
753 
754  ! ... Local variables
755  INTEGER(4) :: ND, nAuxVars, jj, iDim, iDir
756  INTEGER(8) :: Nc, l0, iPoint, pointIndex, metricOffset
757  INTEGER(4) :: normDir, sgn ! , gas_dv_model
758  INTEGER(4) :: numDebugPoints, iDebugPoint,debugPoints(6)
759  ! INTEGER, Dimension(MAX_ND) :: N, Np
760  ! INTEGER, Pointer :: iblank(:), t2Map(:,:)
761  ! INTEGER, Pointer :: grid_is(:), grid_ie(:), patch_is(:), patch_ie(:)
762  ! TYPE (t_grid), Pointer :: grid
763  ! TYPE (t_mixt), Pointer :: state
764  ! TYPE (t_mixt_input), Pointer :: input
765  REAL(8) :: dsgn, sndspdref2, invtempref, tempref, gamref, spcGasConst
766  REAL(8) :: XI_X, XI_Y, XI_Z, XI_T, bndry_h,sbpBoundaryWeight, bcWallTemp
767  REAL(8) :: XI_X_TILDE, XI_Y_TILDE, XI_Z_TILDE, gridJacobian
768  REAL(8) :: rhoFac, penaltyFac, scalarPenaltyFac, scalarDiff
769 
770  REAL(8), POINTER :: XIX(:), bVelocity(:), vWall(:), vWallTarget(:), vTarget(:)
771  ! REAL(KIND=RFREAL), Pointer :: XI_TAU(:,:), cv(:,:), gv(:,:)
772  REAL(8), POINTER :: pnorm_vec(:) ! , rhs(:,:)
773  ! REAL(KIND=RFREAL), Pointer :: patch_ViscousFlux(:,:,:)
774  ! REAL(KIND=RFREAL), Pointer :: patch_ViscousFluxAux(:,:,:)
775  ! REAL(KIND=RFREAL), Pointer :: patch_ViscousRHSAux(:,:,:), rhs_auxVars(:,:)
776  ! REAL(KIND=RFREAL), Pointer :: auxVars(:,:), Lambda(:,:), ibfac(:)
777  REAL(8), POINTER :: Lambda(:,:), gI1(:), uB(:), Tinv(:,:), Tmat(:,:)
778  ! REAL(KIND=RFREAL), Pointer :: gI1(:), ub(:), Tinv(:,:), Tmat(:,:)
779  REAL(8), POINTER :: Aprime(:,:), gI2(:), matX(:,:), penalty(:) ,correction(:)
780  ! REAL(KIND=RFREAL), Pointer :: auxVarsTarget(:,:), cvTarget(:,:)
781  ! REAL(KIND=RFREAL), Pointer :: SBP_invPdiag_block1(:), JAC(:), MT1(:,:)
782  REAL(8) :: norm_vec(3), gamma, uCon, ubAux, specRad, vDotXi, xh
783  REAL(8) :: sigmaI1, sigmaI2, sigmaAux, Sfactor, ibfac_local, T_wall
784  REAL(8) :: spd_snd, RE, REInv, SAT_sigmaI1, SAT_sigmaI2, pressureTarget
785  REAL(8) :: viscousPenaltyScale
786 
787  LOGICAL :: TRIGGERED
788  triggered = .false.
789 
790  ! ... Main data containers
791  ! patch_gridID = patch%gridID
792  ! grid => region%grid(patch_gridID)
793  ! state => region%state(patch_gridID)
794  ! input => grid%input
795 
796  ! ... Pointer dereferences
797  ! grid_is => grid%is
798  ! grid_ie => grid%ie
799  ! patch_is => patch%is
800  ! patch_ie => patch%ie
801  ! cv => state%cv
802  ! gv => state%gv
803  ! cvTarget => state%cvTarget
804  ! rhs => state%rhs
805  ! iblank => grid%iblank
806  ! JAC => grid%JAC
807  ! MT1 => grid%MT1
808  ! ibfac => grid%ibfac
809  ! t2Map => region%global%t2Map
810  ! XI_TAU => grid%XI_TAU
811  ! rhs_auxVars => state%rhs_auxVars
812  ! patch_ViscousFlux => patch%ViscousFlux
813  ! patch_ViscousFluxAux => patch%ViscousFluxAux
814  ! SBP_invPdiag_block1 => grid%SBP_invPdiag_block1
815  ! auxVars => state%auxVars
816  ! auxVarsTarget => state%auxVarsTarget
817 
818  ! ... Problem size
819  nd = numdim
820  nc = numpointsbuffer
821 
822  nauxvars = numscalar
823 
824  ! ... Useful quantities about this boundary
825  normdir = abs(patchnormaldir)
826  sgn = normdir / patchnormaldir
827  dsgn = dble(sgn)
828  numdebugpoints = 0
829 
830  ! ... Reference quantities
831  ! gasParams = (sndspdref,gamma,Cp,1/Re)
832  ! bcParams(:) = sigma1 sigma2
833  ! MJA changed for new nonDimensionalization
834  !sndspdref2 = gasParams(1)*gasParams(1) ! input%sndspdref * input%sndspdref
835  !invtempref = gasParams(3)/sndspdref2 ! input%Cpref / sndspdref2
836  !tempref = 1.0_8 / invtempref
837  !gamref = gasParams(2)
838  !REinv = gasParams(4)
839  gamref = gasparams(2)
840  spcgasconst= gasparams(3)
841  sndspdref2 = gasparams(5)*gasparams(5)
842  tempref = gasparams(4)
843  invtempref = tempref
844  re = gasparams(1)
845  reinv = 1.0_8/re
846  viscouspenaltyscale = gasparams(6)
847 
848  ! ... BC Constants
849  ! bcParams = (sigma1,sbpBoundaryStencilWeight)
850  sat_sigmai1 = bcparams(1) ! input%SAT_sigmaI1_FF
851  sat_sigmai2 = bcparams(2) ! /input%SAT_sigmaI2
852  bcwalltemp = bcparams(3) ! input%bcic_wallTemp
853  sbpboundaryweight = bcparams(4) !
854 
855  !write(*,*) "SAT_sigmaI1 SAT_sigmaI2 bcWallTemp sbpBoundaryWeight tempRef"
856  !write(*,*) SAT_sigmaI1, SAT_sigmaI2, bcWallTemp, sbpBoundaryWeight, tempRef
857 
858  ! gas_dv_model = input%gas_dv_model
859 
860  ! ... memory allocation
861  ALLOCATE(xix(4),vwall(nd),bvelocity(nd),vwalltarget(nd),vtarget(nd))
862  ALLOCATE(gi1(nd+2), ub(nd+2), tmat(nd+2,nd+2), tinv(nd+2,nd+2), penalty(nd+2),correction(nd+2))
863  ALLOCATE(gi2(nd+2), aprime(nd+2,nd+2), lambda(nd+2,nd+2), matx(nd+2,nd+2), pnorm_vec(nd))
864 
865  xix(:) = 0.0_8
866  IF(gridtype <= unirect) THEN
867  xix(normdir) = gridmetric(normdir)
868  metricoffset = 0
869  gridjacobian = jacobiandeterminant(1)
870  ELSE IF(gridtype == rectilinear) THEN
871  metricoffset = (normdir-1)*numpointsbuffer
872  ELSE
873  metricoffset = (normdir-1)*numdim*numpointsbuffer
874  ENDIF
875 
876 ! WRITE(*,*) 'XIX = ',XIX
877 ! WRITE(*,*) 'NORMAL DIR: ',patchNormalDir
878 ! WRITE(*,*) 'SNDSPDREF2 = ',sndspdref2
879 ! WRITE(*,*) 'tempref = ',tempref
880 ! WRITE(*,*) 'gamref = ',gamref
881 ! WRITE(*,*) 'sigma1 = ',SAT_sigmaI1
882 ! WRITE(*,*) 'REInv = ',REinv
883 ! WRITE(*,*) 'Bweight = ',sbpBoundaryWeight
884 ! WRITE(*,*) 'JAC = ',jacobianDeterminant
885 
886 ! TRIGGERED = .TRUE.
887  debugpoints(:) = 0
888  numdebugpoints = 6
889  debugpoints(1) = 91427
890  debugpoints(2) = 91428
891  debugpoints(3) = 91678
892  debugpoints(4) = 114268
893  debugpoints(5) = 114269
894  debugpoints(6) = 114018
895 
896  ! ... loop over all of the points
897  DO ipoint = 1, numpatchpointsop
898 
899  l0 = patchpointsop(ipoint) + 1 ! assume the points coming in are from C (0-based)
900  ! lp = (k-patch_is(3))*Np(1)*Np(2) + (j-patch_is(2))*Np(1) + (i-patch_is(1)+1)
901 
902  ! IF(iPoint == (numPatchPointsOp/2)) TRIGGERED = .FALSE.
903 
904  ! ... iblank factor (either 1 or zero)
905  ibfac_local = 1.0_8 ! ibfac(l0)
906 
907 
908  ! ... construct the normalized metrics
909  ! XI_X = MT1(l0,t2Map(normDir,1)); XI_Y = 0.0_8; XI_Z = 0.0_8
910  ! if (ND >= 2) XI_Y = MT1(l0,t2Map(normDir,2))
911  ! if (ND == 3) XI_Z = MT1(l0,t2Map(normDir,3))
912  ! XI_T = XI_TAU(l0,normDir)
913  ! IF(ND > 1) XI_Y = gridMetric(2)
914  ! IF(ND > 2) XI_Z = gridMetric(3)
915  IF(gridtype == curvilinear ) THEN
916  DO idir = 1,numdim
917  xix(idir) = gridmetric(metricoffset + (idir-1)*numpointsbuffer +l0)
918  END DO
919  gridjacobian = jacobiandeterminant(l0)
920  ELSE IF(gridtype == rectilinear) THEN
921  xix(normdir) = gridmetric(metricoffset + l0)
922  gridjacobian = jacobiandeterminant(l0)
923  ENDIF
924 
925 
926  ! ... multiply by the Jacobian
927  xi_x = xix(1) * gridjacobian ! JAC(l0)
928  xi_y = xix(2) * gridjacobian ! JAC(l0)
929  xi_z = xix(3) * gridjacobian ! JAC(l0)
930  xi_t = xix(4) * gridjacobian ! JAC(l0)
931 
932  ! ... normalized metrics
933  xh = sqrt( (xi_x * xi_x) + (xi_y * xi_y) + (xi_z * xi_z) )
934  bndry_h = 1.0_8 / xh
935  xi_x_tilde = xi_x * bndry_h
936  xi_y_tilde = xi_y * bndry_h
937  xi_z_tilde = xi_z * bndry_h
938 
939 
940  ! ... follow Magnus's algorithm (norm_vec has unit length)
941  norm_vec(1:3) = (/ xi_x_tilde, xi_y_tilde, xi_z_tilde /)
942  pnorm_vec(:) = norm_vec(1:nd) * xh
943 
944  ! ... penalty parameter
945  sigmai1 = -sat_sigmai1 * dsgn
946 
947  ! ... inviscid penalty parameter
948  ! sigmaI1 = -SAT_sigmaI1_FF * dsgn
949 
950  ! ... pull off the boundary data
951  ! uB(:) = cv(l0,:)
952  ub(1) = rhobuffer(l0)
953  vdotxi = 0
954  DO idim = 1, numdim
955  ub(1+idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)
956  bvelocity(idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)/rhobuffer(l0)
957 ! vDotXi = vDotXi + bVelocity(iDim)*XIX(iDim)*gridJacobian*bndry_h
958  vdotxi = vdotxi + bvelocity(idim)*norm_vec(idim)
959  END DO
960  ub(numdim+2) = rhoebuffer(l0)
961 
962  vwall(1:numdim) = vdotxi*norm_vec(1:numdim)
963  ! NON MOVING GRIDS!!
964  vwalltarget(1:numdim) = 0.0_8
965  ! = dot_product(XYZ_TAU(l0,1:ND),norm_vec(1:ND)) * norm_vec(1:ND)
966 
967 
968  DO idim = 1,numdim
969  vtarget(idim) = bvelocity(idim) - ( vwall(idim) - vwalltarget(idim) )
970  END DO
971 
972 ! pressureTarget = pressureBuffer(l0)
973 
974  ! ... target data
975  ! gI1(:) = cvTarget(l0,:)
976  gi1(1) = ub(1)
977  gi1(numdim+2) = ub(numdim+2)
978  DO idim = 1, numdim
979  gi1(1+idim) = gi1(1) * vtarget(idim)
980  gi1(numdim+2) = gi1(numdim+2) + 0.5_8*gi1(1)*(vtarget(idim)**2 - bvelocity(idim)**2)
981  END DO
982 
983 
984  ! ... use target data to make implicit work better?
985  gamma = gamref ! gv(l0,1)
986  ! If (gas_dv_model == DV_MODEL_IDEALGAS .or. &
987  ! gas_dv_model == DV_MODEL_IDEALGAS_MIXTURE) then
988 
989 
990 ! IF(TRIGGERED .EQV. .FALSE.) THEN
991 ! WRITE(*,*) 'L0 = ',l0
992 ! WRITE(*,*) 'XIX = ',XIX
993 ! WRITE(*,*) '1/xh = ',1.0_8/xh
994 ! WRITE(*,*) 'VdotXi = ',vdotxi
995 ! WRITE(*,*) 'bVelocity = ',bVelocity
996 ! WRITE(*,*) 'pnorm: ',pnorm_vec
997 ! WRITE(*,*) 'State: ',uB
998 ! WRITE(*,*) 'Target: ',gI1
999 ! TRIGGERED = .TRUE.
1000 ! ENDIF
1001 
1002 ! DO iDebugPoint = 1, numDebugPoints
1003 ! IF(l0 == debugPoints(iDebugPoint)) THEN
1004 ! WRITE(*,*) 'Before Inviscid point ',l0
1005 ! WRITE(*,*) 'rho = ',ub(1), ' rhoT = ',gI1(1)
1006 ! WRITE(*,*) 'rhoU = ',ub(2),' rhoUT = ',gI1(2)
1007 ! WRITE(*,*) 'rhoV = ',ub(3),' rhoVT = ',gI1(3)
1008 ! WRITE(*,*) 'rhoE = ',ub(4),' rhoET = ',gI1(4)
1009 ! WRITE(*,*) 'U = ',bVelocity(1), ' uT = ',vTarget(1)
1010 ! WRITE(*,*) 'V = ',bVelocity(2), ' uT = ',vTarget(2)
1011 ! WRITE(*,*) 'uWall = ',vWall(1),' uWallTarget = ',vWallTarget(1)
1012 ! WRITE(*,*) 'vWall = ',vWall(2),' vWallTarget = ',vWallTarget(2)
1013 ! WRITE(*,*) 'PNORM = (',pnorm_vec,') vDotXi = ',vDotXi
1014 ! WRITE(*,*) 'normvec = ',norm_vec
1015 ! ENDIF
1016 ! ENDDO
1017 
1018  CALL sat_form_roe_matrices(nd, ub, gi1, gamma, pnorm_vec, tmat, tinv, lambda)
1019 
1020  specrad = 0.0_8
1021  ucon = vdotxi*xh
1022 
1023  ! ... only pick those Lambda that are incoming
1024  If ( sgn == 1 ) Then
1025  Do jj = 1, nd+2
1026  lambda(jj,jj) = max(lambda(jj,jj),0.0_8)
1027  specrad = max(abs(lambda(jj,jj)),abs(specrad))
1028  End Do
1029  Else
1030  Do jj = 1, nd+2
1031  lambda(jj,jj) = min(lambda(jj,jj),0.0_8)
1032  specrad = max(abs(lambda(jj,jj)),abs(specrad))
1033  End Do
1034  End If
1035 
1036 ! IF(TRIGGERED.eqv..FALSE.) THEN
1037 ! WRITE(*,*) 'XI_X_TILDE: ',XI_X_TILDE
1038 ! WRITE(*,*) 'XI_Y_TILDE: ',XI_Y_TILDE
1039 ! WRITE(*,*) 'XI_Z_TILDE: ',XI_Z_TILDE
1040 ! WRITE(*,*) 'XI_X: ',XI_X
1041 ! WRITE(*,*) 'XI_Y: ',XI_Y
1042 ! WRITE(*,*) 'XI_Z: ',XI_Z
1043 ! WRITE(*,*) 'BOUNDARYH: ',bndry_h
1044 ! WRITE(*,*) 'STATE: ',uB(:)
1045 ! WRITE(*,*) 'TARGET:',gI1(:)
1046 ! WRITE(*,*) 'UCON: ',uCon
1047 ! WRITE(*,*) 'SpecRad: ',specRad
1048 ! WRITE(*,*) 'vDotXi: ',vDotXi
1049 ! ENDIF
1050 ! debugPoint = 114273
1051 
1052  matx = matmul(lambda,tinv)
1053 
1054  ! ... compute the characteristic matrix
1055  aprime = matmul(tmat,matx)
1056 
1057  ! ... subtract off the target
1058  ub(1:nd+2) = ub(1:nd+2) - gi1(1:nd+2)
1059 
1060  ! ... compute the characteristic penalty vector
1061  gi2(1:nd+2) = matmul(aprime,ub)
1062 
1063  penaltyfac = sbpboundaryweight*sigmai1
1064  correction(:) = penaltyfac*gi2(:)
1065 
1066 ! DO iDebugPoint = 1, numDebugPoints
1067 ! IF(l0 == debugPoints(iDebugPoint)) THEN
1068 ! WRITE(*,*) 'Before Inviscid point ',l0
1069 ! WRITE(*,*) 'SIGMAI1 = ',sigmaI1,'bWeight = ',sbpBoundaryWeight
1070 ! WRITE(*,*) 'rhoRHS = ',rhoRHS(l0), ' dRho = ',correction(1)
1071 ! WRITE(*,*) 'rhoURHS = ',rhoVRHS(l0),' dRhoU = ',correction(2)
1072 ! WRITE(*,*) 'rhoVRHS = ',rhoVRHS(l0+numPointsBuffer),' dRhoV = ',correction(3)
1073 ! WRITE(*,*) 'rhoERHS = ',rhoERHS(l0),' dRhoE = ',correction(4)
1074 ! WRITE(*,*) 'cpRho = ',gI2(1)
1075 ! WRITE(*,*) 'cpRhoU = ',gI2(2)
1076 ! WRITE(*,*) 'cpRhoV = ',gI2(3)
1077 ! WRITE(*,*) 'cpRhoE = ',gI2(4)
1078 ! ENDIF
1079 ! ENDDO
1080 
1081  ! IBFAC would be considered here
1082  ! Inviscid correction
1083  rhorhs(l0) = rhorhs(l0) + ibfac_local * correction(1)
1084  DO idim = 1,numdim
1085  pointindex = l0 + (idim-1)*numpointsbuffer
1086  rhovrhs(pointindex) = rhovrhs(pointindex) + ibfac_local*correction(1+idim)
1087  END DO
1088  rhoerhs(l0) = rhoerhs(l0) + ibfac_local * correction(numdim+2)
1089 
1090 
1091 ! IF(TRIGGERED.eqv..FALSE.) THEN
1092 ! WRITE(*,*) 'CORRECTION(',l0,'): ',correction(:)
1093 ! TRIGGERED = .TRUE.
1094 ! ENDIF
1095  ! ... Species treatment (inviscid part)
1096  IF(numscalar > 0) THEN
1097  If (dsgn * ucon > 0) Then
1098  sfactor = -ucon * dsgn * ibfac_local * sbpboundaryweight
1099  Do jj = 1, numscalar
1100  pointindex = (jj-1)*numpointsbuffer + l0
1101  scalarrhs(pointindex) = scalarrhs(pointindex) + sfactor*scalarbuffer(pointindex)
1102  End do
1103  End If
1104  END IF
1105 
1106 
1107  ! VISCOUS PART
1108 
1109 ! uB(:) = cv(l0,:)
1110  ! MJA, repeated from above
1111  ub(1) = rhobuffer(l0)
1112  DO idim = 1, numdim
1113  ub(1+idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)
1114  END DO
1115  ub(numdim+2) = rhoebuffer(l0)
1116 
1117  ! ... wall temperature
1118  ! T_wall = bcic_WallTemp / ((gamref - 1.0_8)*tempref)
1119  ! new nondimensionalization...can be moved outside this loop altogether
1120  t_wall = bcwalltemp / tempref
1121 
1122  ! ... target state
1123  ! gI2(1) = RHO_TARGET
1124  gi2(1) = rhobuffer(l0)
1125 
1126  ! ... wall velocity
1127  ! gI2(2:ND+1) = XYZ_TAU(l0,1:ND) * RHO_TARGET
1128  ! gI2(2:ND+1) = 0.0_8
1129  ! MJA, non-moving meshes
1130  gi2(2:idim+1) = 0.0_8
1131 
1132  ! ... thermally perfect or not
1133  ! if (gas_dv_model == DV_MODEL_IDEALGAS) then
1134  ! gI2(ND+2) = RHO_TARGET * T_wall / GAMMA + 0.5_rfreal * SUM(gI2(2:ND+1)**2) / RHO_TARGET
1135 
1136  ! MJA, only implemented gas model right now
1137  ! old nondimen
1138  !gI2(numDim+2) = rhoTarget(l0) * T_wall / gamref + 0.5_8 * SUM(gI2(2:iDim+1)**2) / gI2(1)
1139  ! new dimen, roll spcGasConst and (gamma-1) into one const outside of loop
1140 
1141  ! use a temperature from the input file, if that temperature is <0,
1142  !use the target energy (temperature)
1143  IF (t_wall < 0) THEN
1144  gi2(numdim+2) = rhoetarget(l0)
1145  ELSE
1146  gi2(numdim+2) = gi2(1) * t_wall * spcgasconst / (gamref-1) &
1147  + 0.5_8 * sum(gi2(2:idim+1)**2) / gi2(1)
1148  ENDIF
1149 
1150  ! sigmaI2 = -(SAT_sigmaI2 / 4.0_rfreal) * (SBP_invPdiag_block1(1)/bndry_h)**2 * &
1151  ! maxval(tv(l0,1:)) / RHO_TARGET * MAX(GAMMA / Pr, 5.0_rfreal / 3.0_rfreal)
1152  ! sigmaI2 = sigmaI2 * REinv
1153  sigmai2 = -sat_sigmai2
1154 ! penaltyFac = (sbpBoundaryWeight/bndry_h)**2*sigmaI2/4.0_8 * &
1155 ! REInv*5.0_8/3.0_8*max(muBuffer(l0),lambdaBuffer(l0))/rhoTarget(l0)
1156 ! TESTING MTC
1157  penaltyfac = (sbpboundaryweight/bndry_h)**2*sigmai2/4.0_8 * &
1158  reinv*viscouspenaltyscale*max(mubuffer(l0),lambdabuffer(l0))/gi2(1)
1159 
1160 
1161  ! ... compute penalty
1162  ! ... compute penalty
1163  ! penalty(:) = sigmaI2 * (uB(1:ND+2) - gI2(1:ND+2))
1164  penalty(:) = penaltyfac * (ub(1:nd+2) - gi2(1:nd+2))
1165 
1166  ! ... add to the rhs
1167  ! do jj = 1, ND+2
1168  ! rhs(l0,jj) = rhs(l0,jj) + ibfac_local * penalty(jj)
1169  ! end do
1170  rhorhs(l0) = rhorhs(l0) + ibfac_local * penalty(1)
1171  DO idim = 1,numdim
1172  pointindex = l0 + (idim-1)*numpointsbuffer
1173  rhovrhs(pointindex) = rhovrhs(pointindex) + ibfac_local*penalty(1+idim)
1174  END DO
1175  rhoerhs(l0) = rhoerhs(l0) + ibfac_local * penalty(numdim+2)
1176 
1177 ! DO iDebugPoint = 1, numDebugPoints
1178 ! IF(l0 == debugPoints(iDebugPoint)) THEN
1179 ! WRITE(*,*) 'After Viscous point ',l0
1180 ! WRITE(*,*) 'rhoRHS = ',rhoRHS(l0), ' dRho = ',penalty(1)
1181 ! WRITE(*,*) 'rhoURHS = ',rhoVRHS(l0),' dRhoU = ',penalty(2)
1182 ! WRITE(*,*) 'rhoVRHS = ',rhoVRHS(l0+numPointsBuffer),' dRhoV = ',penalty(3)
1183 ! WRITE(*,*) 'rhoERHS = ',rhoERHS(l0),' dRhoE = ',penalty(4)
1184 ! WRITE(*,*) 'gI2(1) = ',gI2(1)
1185 ! WRITE(*,*) 'gI2(2) = ',gI2(2)
1186 ! WRITE(*,*) 'gI2(3) = ',gI2(3)
1187 ! WRITE(*,*) 'gI2(4) = ',gI2(4)
1188 ! ENDIF
1189 ! ENDDO
1190 
1191  !zero grad for species
1192  ! sigmaAux = ibfac_local * dsgn *1_8* SBP_invPdiag_block1(1) * JAC(l0)
1193  ! do jj = 1, nAuxVars
1194  ! ubAux = patch_ViscousFluxAux(jj,lp,normDir)
1195  ! rhs_auxVars(l0,jj) = rhs_auxVars(l0,jj) + sigmaAux* ubAux
1196  ! end do
1197 
1198  ! MJA
1199  ! species transport not fully implemented yet for viscous terms
1200 
1201  END DO
1202 
1203  DEALLOCATE(xix,vwall,bvelocity,vwalltarget,vtarget)
1204  DEALLOCATE(gi1, ub, tmat, tinv, penalty, correction)
1205  DEALLOCATE(gi2, aprime, lambda, matx, pnorm_vec)
1206 
1207  RETURN
1208 
1209  END SUBROUTINE noslip_isothermal
1210 
1211  SUBROUTINE roe_mat(numDim,inCV,outCV, &
1212  gamma,normDir,tMat,tInv,lambda,status)
1213 
1214  IMPLICIT NONE
1215 
1216  INTEGER(4), INTENT(IN) :: numDim
1217  REAL(8), INTENT(IN) :: inCV(numdim+2), outCV(numdim+2)
1218  REAL(8), INTENT(IN) :: gamma, normDir(numdim)
1219  REAL(8), INTENT(OUT) :: tMat((numdim+2)*(numdim+2)),tInv((numdim+2)*(numdim+2)),lambda(numdim+2)
1220  INTEGER(4), INTENT(OUT) :: status
1221 
1222  REAL(8) :: gammaM1,uCon,alpha,beta,phi2,theta,denomRoe,xiX,xiY,xiZ
1223  REAL(8) :: xiXTilde,xiYTilde,xiZTilde
1224  REAL(8) :: rhoIn,rhoV1In,rhoV2In,rhoV3In,rhoM1In,rhoEIn
1225  REAL(8) :: rhoOut,rhoV1Out,rhoV2Out,rhoV3Out,rhoM1Out,rhoEOut
1226  REAL(8) :: inV1,inV2,inV3,outV1,outV2,outV3
1227  REAL(8) :: inKE,inIE,inP,inH
1228  REAL(8) :: outKE,outIE,outP,outH
1229  REAL(8) :: ccRoe,bbRoe,v1Roe,v2Roe,v3Roe,hRoe,keRoe
1230  REAL(8) :: soundSpeedRoe
1231  INTEGER(4) :: iDim
1232 
1233  status = 0
1234 
1235  IF(numdim == 1 .OR. numdim > 3) THEN
1236  status = 1
1237  RETURN
1238  ENDIF
1239 
1240  ! ... constants
1241  gammam1 = gamma - 1.0_8
1242 
1243  ! form primitives from "in"
1244  rhoin = incv(1)
1245  rhov1in = incv(2)
1246  rhov2in = incv(3)
1247  rhov3in = 0.0_8
1248  rhom1in = 1.0/rhoin
1249 
1250  inv1 = rhom1in*rhov1in
1251  inv2 = rhom1in*rhov2in
1252 
1253  IF(numdim == 2) THEN
1254  rhoein = incv(4)
1255  ELSE
1256  rhov3in = incv(4)
1257  rhoein = incv(5)
1258  ENDIF
1259  inv3 = rhom1in*rhov3in
1260 
1261  inke = 0.5_8 * (inv1**2 + inv2**2 + inv3**2)
1262  inie = rhom1in*rhoein - inke
1263  inp = gammam1*rhoin*inie
1264  inh = (gamma/gammam1)*inp/rhoin + inke
1265 
1266 
1267  ! form primitives from "out"
1268  rhoout = outcv(1)
1269  rhov1out = outcv(2)
1270  rhov2out = outcv(3)
1271  rhom1out = 1.0/rhoout
1272  outv1 = rhom1out*rhov1out
1273  outv2 = rhom1out*rhov2out
1274  IF(numdim == 2) THEN
1275  rhoeout = outcv(4)
1276  ELSE
1277  rhoeout = outcv(5)
1278  rhov3out = outcv(4)
1279  ENDIF
1280  outv3 = rhom1out*rhov3out
1281 
1282  outke = 0.5_8 * (outv1**2 + outv2**2 + outv3**2)
1283  outie = rhom1out*rhoeout - outke
1284  outp = gammam1*rhoout*outie
1285  outh = (gamma/gammam1)*outp/rhoout + outke
1286 
1287  ! Roe vars
1288  ccroe = sqrt(rhoout/rhoin)
1289  bbroe = 1.0_8 / (1.0_8 + ccroe)
1290  v1roe = (inv1 + outv1 * ccroe) * bbroe
1291  v2roe = (inv2 + outv2 * ccroe) * bbroe
1292  v3roe = (inv3 + outv3 * ccroe) * bbroe
1293  hroe = (inh + outh*ccroe)*bbroe
1294  keroe = 0.5_8 * (v1roe**2 + v2roe**2 + v3roe**2)
1295 
1296  IF((hroe - keroe) <= 0.0_8) THEN
1297  write (*,'(5(a,E20.8,1X))') "gamma = ", gamma, ", inKE = ", inke, ", inIE = ", inie, &
1298  ", inP = ", inp, ", inH = ", inh
1299  status = 1
1300  RETURN
1301  ENDIF
1302 
1303  soundspeedroe = sqrt(gammam1*(hroe - keroe))
1304 
1305  ! ... variable renaming
1306  xix = normdir(1)
1307  xiy = normdir(2)
1308  xiz = 0.0_8
1309  IF(numdim == 3) xiz = normdir(3)
1310  denomroe = 1.0_8 / sqrt(xix**2 + xiy**2 + xiz**2)
1311  xixtilde = xix * denomroe
1312  xiytilde = xiy * denomroe
1313  xiztilde = xiz * denomroe
1314 
1315  ucon = v1roe*xix + v2roe*xiy + v3roe*xiz
1316 
1317  ! ... directly from Pulliam & Chaussee
1318  alpha = 0.5_8*rhoin/soundspeedroe
1319  beta = 1.0_8/(rhoin*soundspeedroe)
1320  theta = xixtilde*v1roe + xiytilde*v2roe + xiztilde*v3roe
1321  phi2 = 0.5_8*(gammam1)*(v1roe**2 + v2roe**2 + v3roe**2)
1322 
1323  ! ... compute the diagonal matrix lambda
1324  lambda(:) = 0.0_8
1325  DO idim = 1, numdim
1326  lambda(idim) = ucon
1327  END DO
1328  lambda(numdim+1) = ucon + soundspeedroe * sqrt(xix**2 + xiy**2 + xiz**2)
1329  lambda(numdim+2) = ucon - soundspeedroe * sqrt(xix**2 + xiy**2 + xiz**2)
1330 
1331  if (numdim == 2) then
1332 
1333  tmat(1) = 1.0_8
1334  tmat(2) = v1roe
1335  tmat(3) = v2roe
1336  tmat(4) = phi2/gammam1
1337 
1338  tmat(5) = 0.0_8
1339  tmat(6) = xiytilde*rhoin
1340  tmat(7) = -xixtilde*rhoin
1341  tmat(8) = rhoin*(xiytilde*v1roe - xixtilde*v2roe)
1342 
1343  tmat(9) = alpha
1344  tmat(10) = alpha*(v1roe + xixtilde*soundspeedroe)
1345  tmat(11) = alpha*(v2roe + xiytilde*soundspeedroe)
1346  tmat(12) = alpha*((phi2 + soundspeedroe**2)/gammam1 + soundspeedroe*theta)
1347 
1348  tmat(13) = alpha
1349  tmat(14) = alpha*(v1roe - xixtilde*soundspeedroe)
1350  tmat(15) = alpha*(v2roe - xiytilde*soundspeedroe)
1351  tmat(16) = alpha*((phi2 + soundspeedroe**2)/gammam1 - soundspeedroe*theta)
1352 
1353  tinv(1) = 1.0_8 - phi2 / soundspeedroe**2
1354  tinv(2) = -(xiytilde * v1roe - xixtilde * v2roe) / rhoin
1355  tinv(3) = beta * (phi2 - soundspeedroe * theta)
1356  tinv(4) = beta * (phi2 + soundspeedroe * theta)
1357 
1358  tinv(5) = gammam1 * v1roe / soundspeedroe**2
1359  tinv(6) = xiytilde / rhoin
1360  tinv(7) = beta * (xixtilde * soundspeedroe - gammam1 * v1roe)
1361  tinv(8) = -beta * (xixtilde * soundspeedroe + gammam1 * v1roe)
1362 
1363  tinv(9) = gammam1 * v2roe / soundspeedroe**2
1364  tinv(10) = -xixtilde / rhoin
1365  tinv(11) = beta * (xiytilde * soundspeedroe - gammam1 * v2roe)
1366  tinv(12) = -beta * (xiytilde * soundspeedroe + gammam1 * v2roe)
1367 
1368  tinv(13) = -gammam1 / soundspeedroe**2
1369  tinv(14) = 0.0_8
1370  tinv(15) = beta * gammam1
1371  tinv(16) = beta * gammam1
1372 
1373  ELSE IF (numdim == 3) THEN
1374 
1375 
1376  tmat(1) = xixtilde
1377  tmat(2) = xixtilde * v1roe
1378  tmat(3) = xixtilde * v2roe + xiztilde * rhoin
1379  tmat(4) = xixtilde * v3roe - xiytilde * rhoin
1380  tmat(5) = xixtilde * phi2 / gammam1 + rhoin * (xiztilde * v2roe - xiytilde * v3roe)
1381 
1382  tmat(6) = xiytilde
1383  tmat(7) = xiytilde * v1roe - xiztilde * rhoin
1384  tmat(8) = xiytilde * v2roe
1385  tmat(9) = xiytilde * v3roe + xixtilde * rhoin
1386  tmat(10) = xiytilde * phi2 / gammam1 + rhoin * (xixtilde * v3roe - xiztilde * v1roe)
1387 
1388  tmat(11) = xiztilde
1389  tmat(12) = xiztilde * v1roe + xiytilde * rhoin
1390  tmat(13) = xiztilde * v2roe - xixtilde * rhoin
1391  tmat(14) = xiztilde * v3roe
1392  tmat(15) = xiztilde * phi2 / gammam1 + rhoin * (xiytilde * v1roe - xixtilde * v2roe)
1393 
1394  tmat(16) = alpha
1395  tmat(17) = alpha * (v1roe + xixtilde * soundspeedroe)
1396  tmat(18) = alpha * (v2roe + xiytilde * soundspeedroe)
1397  tmat(19) = alpha * (v3roe + xiztilde * soundspeedroe)
1398  tmat(20) = alpha * ((phi2 + soundspeedroe**2)/gammam1 + soundspeedroe * theta)
1399 
1400  tmat(21) = alpha
1401  tmat(22) = alpha * (v1roe - xixtilde * soundspeedroe)
1402  tmat(23) = alpha * (v2roe - xiytilde * soundspeedroe)
1403  tmat(24) = alpha * (v3roe - xiztilde * soundspeedroe)
1404  tmat(25) = alpha * ((phi2 + soundspeedroe**2)/gammam1 - soundspeedroe * theta)
1405 
1406  tinv(1) = xixtilde * (1.0_8 - phi2 / soundspeedroe**2) - (xiztilde * v2roe - xiytilde * v3roe) / rhoin
1407  tinv(2) = xiytilde * (1.0_8 - phi2 / soundspeedroe**2) - (xixtilde * v3roe - xiztilde * v1roe) / rhoin
1408  tinv(3) = xiztilde * (1.0_8 - phi2 / soundspeedroe**2) - (xiytilde * v1roe - xixtilde * v2roe) / rhoin
1409  tinv(4) = beta * (phi2 - soundspeedroe * theta)
1410  tinv(5) = beta * (phi2 + soundspeedroe * theta)
1411 
1412  tinv(6) = xixtilde * gammam1 * v1roe / soundspeedroe**2
1413  tinv(7) = -xiztilde / rhoin + xiytilde * gammam1 * v1roe / soundspeedroe**2
1414  tinv(8) = xiytilde / rhoin + xiztilde * gammam1 * v1roe / soundspeedroe**2
1415  tinv(9) = beta * (xixtilde * soundspeedroe - gammam1 * v1roe)
1416  tinv(10) = -beta * (xixtilde * soundspeedroe + gammam1 * v1roe)
1417 
1418  tinv(11) = xiztilde / rhoin + xixtilde * gammam1 * v2roe / soundspeedroe**2
1419  tinv(12) = xiytilde * gammam1 * v2roe / soundspeedroe**2
1420  tinv(13) = -xixtilde / rhoin + xiztilde * gammam1 * v2roe / soundspeedroe**2
1421  tinv(14) = beta * (xiytilde * soundspeedroe - gammam1 * v2roe)
1422  tinv(15) = -beta * (xiytilde * soundspeedroe + gammam1 * v2roe)
1423 
1424  tinv(16) = -xiytilde / rhoin + xixtilde * gammam1 * v3roe / soundspeedroe**2
1425  tinv(17) = xixtilde / rhoin + xiytilde * gammam1 * v3roe / soundspeedroe**2
1426  tinv(18) = xiztilde * gammam1 * v3roe / soundspeedroe**2
1427  tinv(19) = beta * (xiztilde * soundspeedroe - gammam1 * v3roe)
1428  tinv(20) = -beta * (xiztilde * soundspeedroe + gammam1 * v3roe)
1429 
1430  tinv(21) = -xixtilde * gammam1 / soundspeedroe**2
1431  tinv(22) = -xiytilde * gammam1 / soundspeedroe**2
1432  tinv(23) = -xiztilde * gammam1 / soundspeedroe**2
1433  tinv(24) = beta * gammam1
1434  tinv(25) = beta * gammam1
1435 
1436 
1437  END IF
1438 
1439  END SUBROUTINE roe_mat
1440 
1441 
1442 
1443 
1444  subroutine sat_form_roe_matrices_tp(u_in, u_out, ND, gamma, gamref, norm, tmat, tinv, lambda)
1445 
1446 ! USE ModGlobal
1447 ! USE ModMPI
1448  IMPLICIT NONE
1449 
1450  Real(8), Pointer :: u_in(:), u_out(:), tmat(:,:), tinv(:,:), norm(:), lambda(:,:)
1451  Real(8) :: gamma, rho_in, ux_in, uy_in, uz_in, ke_in, p_in, en_in, h_in
1452  Real(8) :: gm1, rho_out, ux_out, uy_out, uz_out, ke_out, p_out, en_out, h_out
1453  Real(8) :: cc, bb, ux_Roe, uy_Roe, uz_Roe, h_Roe, ke_Roe, spd_snd_Roe, T_Roe, en_Roe
1454  Real(8) :: XI_X_TILDE, XI_Y_TILDE, XI_Z_TILDE, ucon, rho, alpha, theta, beta, phi2
1455  Real(8) :: XI_X, XI_Y, XI_Z, XI_T, denom, GAMREF, Up, Um
1456 
1457  Integer :: N, ND, jj, i, j, counter,ii
1458  Logical, dimension(ND+2) :: BWORK
1459  Integer, dimension(ND+2) :: order
1460 
1461  ! ... Lapack arguments
1462  CHARACTER(len=1) :: JOBVL, JOBVR
1463  INTEGER :: INFO, LDA, LDVL, LDVR, LWORK,M
1464  Real(8), dimension(ND+2,ND+2) :: A, VL, VR,B,LAM, Pmat
1465  Real(8), dimension(ND+2) :: WI, WR
1466  Real(8), dimension(300) :: WORK
1467  ! ... for right vector inversion
1468  integer, dimension(ND+2) :: IPIV
1469 
1470  gamma = gamref
1471 
1472  ! ... temporary
1473  xi_t = 0.0_8
1474 
1475  jobvl = 'N'
1476  jobvr = 'V'
1477  lda = nd + 2
1478  ldvl = nd + 2
1479  ldvr = nd + 2
1480  m = nd+2
1481 
1482  lwork = 300
1483 
1484  lam(:,:) = 0.0_8
1485  a(:,:) = 0.0_8
1486  vl(:,:) = 0.0_8
1487  vr(:,:) = 0.0_8
1488  wr(:) = 0.0_8
1489  wi(:) = 0.0_8
1490  info = 0
1491 
1492 
1493  ! ... constants
1494  gm1 = gamma - 1.0_8
1495 
1496  ! ... form primative variables from "in"
1497  rho_in = u_in(1)
1498  ux_in = u_in(2) / rho_in; uy_in = 0.0_8; uz_in = 0.0_8
1499  if (nd >= 2) uy_in = u_in(3) / rho_in
1500  if (nd == 3) uz_in = u_in(4) / rho_in
1501  ke_in = 0.5_8 * (ux_in**2 + uy_in**2 + uz_in**2)
1502  en_in = u_in(nd+2)/rho_in - ke_in
1503  p_in = gm1 * rho_in * en_in
1504  h_in = (gamma / gm1) * p_in / rho_in + ke_in
1505 
1506  ! ... form primative variables from "out"
1507  rho_out = u_out(1)
1508  ux_out = u_out(2) / rho_out; uy_out = 0.0_8; uz_out = 0.0_8
1509  if (nd >= 2) uy_out = u_out(3) / rho_out
1510  if (nd == 3) uz_out = u_out(4) / rho_out
1511  ke_out = 0.5_8 * (ux_out**2 + uy_out**2 + uz_out**2)
1512  en_out = u_out(nd+2)/rho_out - ke_out
1513  p_out = gm1 * rho_out * en_out
1514  h_out = (gamma / gm1) * p_out / rho_out + ke_out
1515 
1516  ! ... form Roe variables
1517  cc = sqrt(rho_out / rho_in)
1518  bb = 1.0_8 / (1.0_8 + cc)
1519  ux_roe = (ux_in + ux_out * cc) * bb
1520  uy_roe = (uy_in + uy_out * cc) * bb
1521  uz_roe = (uz_in + uz_out * cc) * bb
1522  h_roe = ( h_in + h_out * cc) * bb
1523  ke_roe = 0.5_8 * (ux_roe**2 + uy_roe**2 + uz_roe**2)
1524  spd_snd_roe = sqrt(gamma/gamref * gm1*(h_roe - ke_roe))
1525 
1526  ! ... temporarily
1527  t_roe = gamma*(en_in + en_out * cc) * bb
1528  en_roe = t_roe/gamma
1529 
1530  ! ... variable renaming
1531  xi_x = norm(1); xi_y = 0.0_8; xi_z = 0.0_8
1532  if (nd >= 2) xi_y = norm(2)
1533  if (nd == 3) xi_z = norm(3)
1534  denom = 1.0_8 / sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1535  xi_x_tilde = xi_x * denom
1536  xi_y_tilde = xi_y * denom
1537  xi_z_tilde = xi_z * denom
1538  rho = rho_in
1539 
1540  ! ... directly from Pulliam & Chaussee
1541  theta = xi_x * ux_roe + xi_y * uy_roe + xi_z * uz_roe
1542  phi2 = 0.5_8 * (gamma - 1.0_8) * (ux_roe**2 + uy_roe**2 + uz_roe**2)
1543 
1544  if (nd == 2) then
1545 
1546  a(1,1) = xi_t
1547  a(1,2) = xi_x
1548  a(1,3) = xi_y
1549  a(1,4) = 0.0_8
1550 
1551  a(2,1) = xi_x * ( t_roe / gamref * (gamref - 1.0_8) - (gm1 * en_roe - phi2)) - theta * ux_roe
1552  a(2,2) = xi_t + theta - xi_x * (gamma - 2.0_8) * ux_roe
1553  a(2,3) = xi_y * ux_roe - gm1 * xi_x * uy_roe
1554  a(2,4) = gm1 * xi_x
1555 
1556  a(3,1) = xi_y * ( t_roe / gamref * (gamref - 1.0_8) - (gm1 * en_roe - phi2)) - theta * uy_roe
1557  a(3,2) = xi_x * uy_roe - gm1 * xi_y * ux_roe
1558  a(3,3) = xi_t + theta - xi_y * (gamma - 2.0_8) * uy_roe
1559  a(3,4) = gm1 * xi_y
1560 
1561  a(4,1) = theta * (2.0_8 * phi2 - gamma * (en_roe + ke_roe))
1562  a(4,2) = xi_x * (t_roe / gamref * (gamref - 1.0_8) + en_roe + phi2 / gm1 ) - gm1 * theta * ux_roe
1563  a(4,3) = xi_y * (t_roe / gamref * (gamref - 1.0_8) + en_roe + phi2 / gm1 ) - gm1 * theta * uy_roe
1564  a(4,4) = xi_t + gamma * theta
1565 
1566  else if (nd == 3) then
1567 
1568  a(1,1) = xi_t
1569  a(1,2) = xi_x
1570  a(1,3) = xi_y
1571  a(1,4) = xi_z
1572  a(1,5) = 0.0_8
1573 
1574  a(2,1) = xi_x * ( t_roe / gamref * (gamref - 1.0_8) - (gm1 * en_roe - phi2)) - theta * ux_roe
1575  a(2,2) = xi_t + theta - xi_x * (gamma - 2.0_8) * ux_roe
1576  a(2,3) = xi_y * ux_roe - gm1 * xi_x * uy_roe
1577  a(2,4) = xi_z * ux_roe - gm1 * xi_x * uz_roe
1578  a(2,5) = gm1 * xi_x
1579 
1580  a(3,1) = xi_y * ( t_roe / gamref * (gamref - 1.0_8) - (gm1 * en_roe - phi2)) - theta * uy_roe
1581  a(3,2) = xi_x * uy_roe - gm1 * xi_y * ux_roe
1582  a(3,3) = xi_t + theta - xi_y * (gamma - 2.0_8) * uy_roe
1583  a(3,4) = xi_z * uy_roe - gm1 * xi_y * uz_roe
1584  a(3,5) = gm1 * xi_y
1585 
1586  a(4,1) = xi_z * ( t_roe / gamref * (gamref - 1.0_8) - (gm1 * en_roe - phi2)) - theta * uz_roe
1587  a(4,2) = xi_x * uz_roe - gm1 * xi_z * ux_roe
1588  a(4,3) = xi_y * uz_roe - gm1 * xi_z * uy_roe
1589  a(4,4) = xi_t + theta - xi_z * (gamma - 2.0_8) * uz_roe
1590  a(4,5) = gm1 * xi_z
1591 
1592  a(5,1) = theta * (2.0_8 * phi2 - gamma * (en_roe + ke_roe))
1593  a(5,2) = xi_x * (t_roe / gamref * (gamref - 1.0_8) + en_roe + phi2 / gm1 ) - gm1 * theta * ux_roe
1594  a(5,3) = xi_y * (t_roe / gamref * (gamref - 1.0_8) + en_roe + phi2 / gm1 ) - gm1 * theta * uy_roe
1595  a(5,4) = xi_z * (t_roe / gamref * (gamref - 1.0_8) + en_roe + phi2 / gm1 ) - gm1 * theta * uz_roe
1596  a(5,5) = xi_t + gamma * theta
1597 
1598  end if
1599 
1600  ! ... find the Eigenvalues and Eigenvectors
1601  ! ... Lapack driver DGEEV
1602 !#ifdef HAVE_BLAS
1603 ! CALL DGEEV( JOBVL, JOBVR, M, A, LDA, WR, WI, VL, LDVL, VR, &
1604 ! LDVR, WORK, LWORK, INFO )
1605 !#endif
1606 
1607  lambda(:,:) = 0.0_8
1608  pmat(:,:) = 0.0_8
1609 
1610  ! ... form transformation matrices
1611  ucon = ux_roe * xi_x + uy_roe * xi_y + uz_roe * xi_z
1612  up = ucon + 0.1_8 * spd_snd_roe * sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1613  um = ucon - 0.1_8 * spd_snd_roe * sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1614  counter = 0
1615 
1616  ! ... Sort eigenvalues: {U,U,U,U+c,U-c}
1617  do jj = 1, nd+2
1618 
1619  if(wr(jj) .lt. up .and. wr(jj) .gt. um) then
1620 
1621  counter = counter + 1
1622  lambda(counter,counter) = wr(jj)
1623  order(jj) = counter
1624  pmat(jj,counter) = 1.0_8
1625  elseif(wr(jj) .gt. ucon) then
1626  lambda(nd+1,nd+1) = wr(jj)
1627  order(jj) = nd+1
1628  pmat(jj,nd+1) = 1.0_8
1629  elseif(wr(jj) .lt. ucon) then
1630  lambda(nd+2,nd+2) = wr(jj)
1631  order(jj) = nd+2
1632  pmat(jj,nd+2) = 1.0_8
1633  end if
1634  end do
1635 
1636  tmat = matmul(vr,pmat)
1637  vr = tmat
1638 
1639  ! ... invert VR
1640  ! ... using a pseudoinvert using singular value decomposition
1641  ! ... VR = U*D*V^T
1642 
1643  ! ... Args: A,A, size VR, size VR, VR, size VR, array of singular values WI, U, size U, V^T, size V, work, lwork, INFO
1644  jobvl = 'A'
1645  wi(:) = 0.0_8
1646  vl(:,:) = 0.0_8
1647  a(:,:) = 0.0_8
1648 
1649  ! ... VR = VL*diag(WI)*A^T
1650 !#ifdef HAVE_BLAS
1651 ! call DGESVD(JOBVL,JOBVL,LDVR,LDVR,VR,LDVR,WI,VL,LDVR,A,LDVR,WORK,LWORK,INFO)
1652 !#endif
1653 
1654  tinv(:,:) = 0.0_8
1655  ! ... pseudo inverse is inv(VR) = inv(VL*diag(WI)*A^T) = A*inv(diag(WI))*VL^T where small (~0) values of WI -> 1/(WI) = 0
1656  do i = 1,nd+2
1657  if(wi(i) .gt. tiny(wi(i))) then
1658  tinv(i,i) = 1.0_8/wi(i)
1659  else
1660  tinv(i,i) = 0.0_8
1661  end if
1662  end do
1663 
1664 
1665  tinv = matmul(transpose(a),tinv)
1666  tinv = matmul(tinv,transpose(vl))
1667 
1668  end subroutine sat_form_roe_matrices_tp
1669 
1670 
1671 
1672  subroutine sat_form_roe_matrices(ND, u_in, u_out, gamma, norm, tmat, tinv, lambda)
1673 
1674 ! USE ModGlobal
1675 
1676  IMPLICIT NONE
1677 
1678  Real(8) :: u_in(nd+2), u_out(nd+2), tmat(nd+2,nd+2), tinv(nd+2,nd+2), norm(nd), lambda(nd+2,nd+2)
1679  Real(8) :: gamma, rho_in, ux_in, uy_in, uz_in, ke_in, p_in, en_in, h_in
1680  Real(8) :: gm1, rho_out, ux_out, uy_out, uz_out, ke_out, p_out, en_out, h_out
1681  Real(8) :: cc, bb, ux_Roe, uy_Roe, uz_Roe, h_Roe, ke_Roe, spd_snd_Roe
1682  Real(8) :: XI_X_TILDE, XI_Y_TILDE, XI_Z_TILDE, ucon, rho, alpha, theta, beta, phi2
1683  Real(8) :: XI_X, XI_Y, XI_Z, denom
1684  Integer :: N, ND, jj
1685 
1686  ! ... size of this problem
1687  !ND = size(u_in) - 2
1688 
1689  ! ... constants
1690  gm1 = gamma - 1.0_8
1691 
1692  ! ... form primative variables from "in"
1693  rho_in = u_in(1)
1694  ux_in = u_in(2) / rho_in; uy_in = 0.0_8; uz_in = 0.0_8
1695  if (nd >= 2) uy_in = u_in(3) / rho_in
1696  if (nd == 3) uz_in = u_in(4) / rho_in
1697  ke_in = 0.5_8 * (ux_in**2 + uy_in**2 + uz_in**2)
1698  en_in = u_in(nd+2)/rho_in - ke_in
1699  p_in = gm1 * rho_in * en_in
1700  h_in = (gamma / gm1) * p_in / rho_in + ke_in
1701 
1702  ! ... form primative variables from "out"
1703  rho_out = u_out(1)
1704  ux_out = u_out(2) / rho_out; uy_out = 0.0_8; uz_out = 0.0_8
1705  if (nd >= 2) uy_out = u_out(3) / rho_out
1706  if (nd == 3) uz_out = u_out(4) / rho_out
1707  ke_out = 0.5_8 * (ux_out**2 + uy_out**2 + uz_out**2)
1708  en_out = u_out(nd+2)/rho_out - ke_out
1709  p_out = gm1 * rho_out * en_out
1710  h_out = (gamma / gm1) * p_out / rho_out + ke_out
1711 
1712  ! ... form Roe variables
1713  cc = sqrt(rho_out / rho_in)
1714  bb = 1.0_8 / (1.0_8 + cc)
1715  ux_roe = (ux_in + ux_out * cc) * bb
1716  uy_roe = (uy_in + uy_out * cc) * bb
1717  uz_roe = (uz_in + uz_out * cc) * bb
1718  h_roe = ( h_in + h_out * cc) * bb
1719  ke_roe = 0.5_8 * (ux_roe**2 + uy_roe**2 + uz_roe**2)
1720  if (h_roe - ke_roe <= 0.0_8) then
1721  ! write (*,'(6(E20.8,1X))') gamma, h_Roe, ke_Roe, ux_roe, uy_roe, uz_roe
1722  write (*,'(5(a,E20.8,1X))') "SATUTIL:FORM_ROE_MATRICES: ERRONEOUS CONDITIONS: gamma = ", gamma, ", ke_in = ", &
1723  ke_in, ", en_in = ", en_in, ", p_in = ", p_in, ", h_in = ", h_in
1724  stop
1725  end if
1726  spd_snd_roe = sqrt(gm1*(h_roe - ke_roe))
1727 
1728  ! ... variable renaming
1729  xi_x = norm(1); xi_y = 0.0_8; xi_z = 0.0_8
1730  if (nd >= 2) xi_y = norm(2)
1731  if (nd == 3) xi_z = norm(3)
1732  denom = 1.0_8 / sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1733  xi_x_tilde = xi_x * denom
1734  xi_y_tilde = xi_y * denom
1735  xi_z_tilde = xi_z * denom
1736  rho = rho_in
1737 
1738 !!$ ! ... debugging
1739 !!$ RHO = rho_in
1740 !!$ UX_ROE = UX_IN
1741 !!$ UY_ROE = UY_IN
1742 !!$ UZ_ROE = UZ_IN
1743 !!$ SPD_SND_ROE = sqrt(gm1*(H_IN - KE_IN))
1744 
1745  ! ... form transformation matrices
1746  ucon = ux_roe * xi_x + uy_roe * xi_y + uz_roe * xi_z
1747 
1748  if (nd == 2) then
1749 
1750  ! ... directly from Pulliam & Chaussee
1751  alpha = 0.5_8 * rho / spd_snd_roe
1752  beta = 1.0_8 / (rho * spd_snd_roe)
1753  theta = xi_x_tilde * ux_roe + xi_y_tilde * uy_roe
1754  phi2 = 0.5_8 * (gamma - 1.0_8) * (ux_roe**2 + uy_roe**2)
1755  tmat(1,1) = 1.0_8
1756  tmat(1,2) = 0.0_8
1757  tmat(1,3) = alpha
1758  tmat(1,4) = alpha
1759  tmat(2,1) = ux_roe
1760  tmat(2,2) = xi_y_tilde * rho
1761  tmat(2,3) = alpha * (ux_roe + xi_x_tilde * spd_snd_roe)
1762  tmat(2,4) = alpha * (ux_roe - xi_x_tilde * spd_snd_roe)
1763  tmat(3,1) = uy_roe
1764  tmat(3,2) = -xi_x_tilde * rho
1765  tmat(3,3) = alpha * (uy_roe + xi_y_tilde * spd_snd_roe)
1766  tmat(3,4) = alpha * (uy_roe - xi_y_tilde * spd_snd_roe)
1767  tmat(4,1) = phi2 / (gamma - 1.0_8)
1768  tmat(4,2) = rho * (xi_y_tilde * ux_roe - xi_x_tilde * uy_roe)
1769  tmat(4,3) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) + spd_snd_roe * theta)
1770  tmat(4,4) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) - spd_snd_roe * theta)
1771 
1772  tinv(1,1) = 1.0_8 - phi2 / spd_snd_roe**2
1773  tinv(1,2) = (gamma - 1.0_8) * ux_roe / spd_snd_roe**2
1774  tinv(1,3) = (gamma - 1.0_8) * uy_roe / spd_snd_roe**2
1775  tinv(1,4) = -(gamma - 1.0_8) / spd_snd_roe**2
1776  tinv(2,1) = -(xi_y_tilde * ux_roe - xi_x_tilde * uy_roe) / rho
1777  tinv(2,2) = xi_y_tilde / rho
1778  tinv(2,3) = -xi_x_tilde / rho
1779  tinv(2,4) = 0.0_8
1780  tinv(3,1) = beta * (phi2 - spd_snd_roe * theta)
1781  tinv(3,2) = beta * (xi_x_tilde * spd_snd_roe - (gamma - 1.0_8) * ux_roe)
1782  tinv(3,3) = beta * (xi_y_tilde * spd_snd_roe - (gamma - 1.0_8) * uy_roe)
1783  tinv(3,4) = beta * (gamma - 1.0_8)
1784  tinv(4,1) = beta * (phi2 + spd_snd_roe * theta)
1785  tinv(4,2) = -beta * (xi_x_tilde * spd_snd_roe + (gamma - 1.0_8) * ux_roe)
1786  tinv(4,3) = -beta * (xi_y_tilde * spd_snd_roe + (gamma - 1.0_8) * uy_roe)
1787  tinv(4,4) = beta * (gamma - 1.0_8)
1788 
1789  ! ... compute the diagonal matrix lambda
1790  lambda(:,:) = 0.0_8
1791  do jj = 1, nd
1792  lambda(jj,jj) = ucon
1793  end do
1794  lambda(nd+1,nd+1) = ucon + spd_snd_roe * sqrt(xi_x**2 + xi_y**2)
1795  lambda(nd+2,nd+2) = ucon - spd_snd_roe * sqrt(xi_x**2 + xi_y**2)
1796 
1797  else if (nd == 3) then
1798 
1799  ! ... directly from Pulliam & Chaussee
1800  alpha = 0.5_8 * rho / spd_snd_roe
1801  beta = 1.0_8 / (rho * spd_snd_roe)
1802  theta = xi_x_tilde * ux_roe + xi_y_tilde * uy_roe + xi_z_tilde * uz_roe
1803  phi2 = 0.5_8 * (gamma - 1.0_8) * (ux_roe**2 + uy_roe**2 + uz_roe**2)
1804 
1805  tmat(1,1) = xi_x_tilde
1806  tmat(1,2) = xi_y_tilde
1807  tmat(1,3) = xi_z_tilde
1808  tmat(1,4) = alpha
1809  tmat(1,5) = alpha
1810  tmat(2,1) = xi_x_tilde * ux_roe
1811  tmat(2,2) = xi_y_tilde * ux_roe - xi_z_tilde * rho
1812  tmat(2,3) = xi_z_tilde * ux_roe + xi_y_tilde * rho
1813  tmat(2,4) = alpha * (ux_roe + xi_x_tilde * spd_snd_roe)
1814  tmat(2,5) = alpha * (ux_roe - xi_x_tilde * spd_snd_roe)
1815  tmat(3,1) = xi_x_tilde * uy_roe + xi_z_tilde * rho
1816  tmat(3,2) = xi_y_tilde * uy_roe
1817  tmat(3,3) = xi_z_tilde * uy_roe - xi_x_tilde * rho
1818  tmat(3,4) = alpha * (uy_roe + xi_y_tilde * spd_snd_roe)
1819  tmat(3,5) = alpha * (uy_roe - xi_y_tilde * spd_snd_roe)
1820  tmat(4,1) = xi_x_tilde * uz_roe - xi_y_tilde * rho
1821  tmat(4,2) = xi_y_tilde * uz_roe + xi_x_tilde * rho
1822  tmat(4,3) = xi_z_tilde * uz_roe
1823  tmat(4,4) = alpha * (uz_roe + xi_z_tilde * spd_snd_roe)
1824  tmat(4,5) = alpha * (uz_roe - xi_z_tilde * spd_snd_roe)
1825  tmat(5,1) = xi_x_tilde * phi2 / (gamma - 1.0_8) + rho * (xi_z_tilde * uy_roe - xi_y_tilde * uz_roe)
1826  tmat(5,2) = xi_y_tilde * phi2 / (gamma - 1.0_8) + rho * (xi_x_tilde * uz_roe - xi_z_tilde * ux_roe)
1827  tmat(5,3) = xi_z_tilde * phi2 / (gamma - 1.0_8) + rho * (xi_y_tilde * ux_roe - xi_x_tilde * uy_roe)
1828  tmat(5,4) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) + spd_snd_roe * theta)
1829  tmat(5,5) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) - spd_snd_roe * theta)
1830 
1831  tinv(1,1) = xi_x_tilde * (1.0_8 - phi2 / spd_snd_roe**2) - (xi_z_tilde * uy_roe - xi_y_tilde * uz_roe) / rho
1832  tinv(1,2) = xi_x_tilde * (gamma - 1.0_8) * ux_roe / spd_snd_roe**2
1833  tinv(1,3) = xi_z_tilde / rho + xi_x_tilde * (gamma - 1.0_8) * uy_roe / spd_snd_roe**2
1834  tinv(1,4) = -xi_y_tilde / rho + xi_x_tilde * (gamma - 1.0_8) * uz_roe / spd_snd_roe**2
1835  tinv(1,5) = -xi_x_tilde * (gamma - 1.0_8) / spd_snd_roe**2
1836  tinv(2,1) = xi_y_tilde * (1.0_8 - phi2 / spd_snd_roe**2) - (xi_x_tilde * uz_roe - xi_z_tilde * ux_roe) / rho
1837  tinv(2,2) = -xi_z_tilde / rho + xi_y_tilde * (gamma - 1.0_8) * ux_roe / spd_snd_roe**2
1838  tinv(2,3) = xi_y_tilde * (gamma - 1.0_8) * uy_roe / spd_snd_roe**2
1839  tinv(2,4) = xi_x_tilde / rho + xi_y_tilde * (gamma - 1.0_8) * uz_roe / spd_snd_roe**2
1840  tinv(2,5) = -xi_y_tilde * (gamma - 1.0_8) / spd_snd_roe**2
1841  tinv(3,1) = xi_z_tilde * (1.0_8 - phi2 / spd_snd_roe**2) - (xi_y_tilde * ux_roe - xi_x_tilde * uy_roe) / rho
1842  tinv(3,2) = xi_y_tilde / rho + xi_z_tilde * (gamma - 1.0_8) * ux_roe / spd_snd_roe**2
1843  tinv(3,3) = -xi_x_tilde / rho + xi_z_tilde * (gamma - 1.0_8) * uy_roe / spd_snd_roe**2
1844  tinv(3,4) = xi_z_tilde * (gamma - 1.0_8) * uz_roe / spd_snd_roe**2
1845  tinv(3,5) = -xi_z_tilde * (gamma - 1.0_8) / spd_snd_roe**2
1846  tinv(4,1) = beta * (phi2 - spd_snd_roe * theta)
1847  tinv(4,2) = beta * (xi_x_tilde * spd_snd_roe - (gamma - 1.0_8) * ux_roe)
1848  tinv(4,3) = beta * (xi_y_tilde * spd_snd_roe - (gamma - 1.0_8) * uy_roe)
1849  tinv(4,4) = beta * (xi_z_tilde * spd_snd_roe - (gamma - 1.0_8) * uz_roe)
1850  tinv(4,5) = beta * (gamma - 1.0_8)
1851  tinv(5,1) = beta * (phi2 + spd_snd_roe * theta)
1852  tinv(5,2) = -beta * (xi_x_tilde * spd_snd_roe + (gamma - 1.0_8) * ux_roe)
1853  tinv(5,3) = -beta * (xi_y_tilde * spd_snd_roe + (gamma - 1.0_8) * uy_roe)
1854  tinv(5,4) = -beta * (xi_z_tilde * spd_snd_roe + (gamma - 1.0_8) * uz_roe)
1855  tinv(5,5) = beta * (gamma - 1.0_8)
1856 
1857  ! ... compute the diagonal matrix lambda
1858  lambda(:,:) = 0.0_8
1859  do jj = 1, nd
1860  lambda(jj,jj) = ucon
1861  end do
1862  lambda(nd+1,nd+1) = ucon + spd_snd_roe * sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1863  lambda(nd+2,nd+2) = ucon - spd_snd_roe * sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1864 
1865  end if
1866 
1867  end subroutine sat_form_roe_matrices
1868 
1869  SUBROUTINE dissipationweight(numDim,bufferSize,numPoints,bufferInterval, &
1870  sigmaDissipation,sigmaDilatation,dilatationRamp,dilatationCutoff, &
1871  divV,sigmaDiss)
1872 
1873  IMPLICIT NONE
1874 
1875  INTEGER(KIND=4), INTENT(IN) :: numDim
1876  INTEGER(KIND=8), INTENT(IN) :: bufferSize(numdim), numPoints
1877  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
1878  REAL(KIND=8), INTENT(IN) :: sigmaDissipation,sigmaDilatation
1879  REAL(KIND=8), INTENT(IN) :: dilatationRamp,dilatationCutoff
1880  REAL(KIND=8), INTENT(IN) :: divV(numpoints)
1881  REAL(KIND=8), INTENT(OUT) :: sigmaDiss(numpoints)
1882 
1883  INTEGER(KIND=8) :: I, J, K
1884  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1885  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1886 
1887  istart = bufferinterval(1)
1888  iend = bufferinterval(2)
1889  xsize = buffersize(1)
1890 
1891  IF(numdim == 1) THEN
1892  DO i = istart, iend
1893  sigmadiss(i) = (sigmadissipation + &
1894  sigmadilatation*0.5_8*(1.0_8 + tanh(dilatationramp*(dilatationcutoff-divv(i)))))
1895  END DO
1896  ELSE IF(numdim == 2) THEN
1897  jstart = bufferinterval(3)
1898  jend = bufferinterval(4)
1899  DO j = jstart, jend
1900  yindex = (j-1)*xsize
1901  DO i = istart, iend
1902  bufferindex = yindex + i
1903  sigmadiss(bufferindex) = (sigmadissipation + &
1904  sigmadilatation*0.5_8*(1.0_8 + tanh(dilatationramp*(dilatationcutoff-divv(bufferindex)))))
1905  END DO
1906  END DO
1907  ELSE IF(numdim == 3) THEN
1908  jstart = bufferinterval(3)
1909  jend = bufferinterval(4)
1910  kstart = bufferinterval(5)
1911  kend = bufferinterval(6)
1912  nplane = xsize * buffersize(2)
1913  DO k = kstart, kend
1914  zindex = (k-1)*nplane
1915  DO j = jstart, jend
1916  yzindex = zindex + (j-1)*xsize
1917  DO i = istart, iend
1918  bufferindex = yzindex + i
1919  sigmadiss(bufferindex) = (sigmadissipation + &
1920  sigmadilatation*0.5_8*(1.0_8 + tanh(dilatationramp*(dilatationcutoff-divv(bufferindex)))))
1921  END DO
1922  END DO
1923  END DO
1924  ENDIF
1925 
1926 
1927  END SUBROUTINE dissipationweight
1928 
1929 END MODULE satutil
satutil::noslip_isothermal
subroutine noslip_isothermal(numDim, bufferSizes, numPointsBuffer, patchNormalDir, patchSizes, numPointsPatch, numPatchPointsOp, patchPointsOp, gridType, gridMetric, jacobianDeterminant, bcParams, gasParams, rhoBuffer, rhoVBuffer, rhoEBuffer, numscalar, scalarBuffer, rhoRHS, rhoVRHS, rhoERHS, scalarRHS, rhoTarget, rhoVTarget, rhoETarget, scalarTarget, muBuffer, lambdaBuffer)
Definition: SATUtil.f90:727
satutil::dissipationweight
subroutine dissipationweight(numDim, bufferSize, numPoints, bufferInterval, sigmaDissipation, sigmaDilatation, dilatationRamp, dilatationCutoff, divV, sigmaDiss)
Definition: SATUtil.f90:1872
satutil::roe_mat
subroutine roe_mat(numDim, inCV, outCV, gamma, normDir, tMat, tInv, lambda, status)
Definition: SATUtil.f90:1213
grid::unirect
integer(kind=4), parameter unirect
Definition: Grid.f90:6
satutil::farfield
subroutine farfield(numDim, bufferSizes, numPointsBuffer, patchNormalDir, patchSizes, numPointsPatch, numPatchPointsOp, patchPointsOp, gridType, gridMetric, jacobianDeterminant, bcParams, gasParams, rhoBuffer, rhoVBuffer, rhoEBuffer, viscousFluxBuffer, numscalar, scalarBuffer, rhoRHS, rhoVRHS, rhoERHS, scalarRHS, rhoTarget, rhoVTarget, rhoETarget, scalarTarget)
Definition: SATUtil.f90:16
grid
Definition: Grid.f90:1
satutil::slip_adiabatic
subroutine slip_adiabatic(numDim, bufferSizes, numPointsBuffer, patchNormalDir, patchSizes, numPointsPatch, numPatchPointsOp, patchPointsOp, gridType, gridMetric, jacobianDeterminant, bcParams, gasParams, rhoBuffer, rhoVBuffer, rhoEBuffer, numscalar, scalarBuffer, rhoRHS, rhoVRHS, rhoERHS, scalarRHS, rhoTarget, rhoVTarget, rhoETarget, scalarTarget)
Definition: SATUtil.f90:360
grid::curvilinear
integer(kind=4), parameter curvilinear
Definition: Grid.f90:8
grid::rectilinear
integer(kind=4), parameter rectilinear
Definition: Grid.f90:7
satutil::sat_form_roe_matrices_tp
subroutine sat_form_roe_matrices_tp(u_in, u_out, ND, gamma, gamref, norm, tmat, tinv, lambda)
Definition: SATUtil.f90:1445
satutil::sat_form_roe_matrices
subroutine sat_form_roe_matrices(ND, u_in, u_out, gamma, norm, tmat, tinv, lambda)
Definition: SATUtil.f90:1673