PlasCom2  1.0
XPACC Multi-physics simluation application
SATUtil.sav.f90
Go to the documentation of this file.
1 MODULE satutil
2 
3  USE grid
4 
5  IMPLICIT NONE
6 
7  TYPE t_nasapolynomials
8  real(8),pointer :: coeffs(:,:)
9  real(8),pointer :: coeffsorig(:,:)
10  real(8) :: hfg
11  real(8),pointer :: cmod(:,:)
12  real(8),pointer :: ch(:,:)
13  real(8),pointer :: cp(:,:)
14  real(8) :: delta
15  end type t_nasapolynomials
16 
17 CONTAINS
18 
19  SUBROUTINE farfield( &
20  numDim,bufferSizes,numPointsBuffer,patchNormalDir,patchSizes, &
21  numPointsPatch,numPatchPointsOp,patchPointsOp,gridType,gridMetric,&
22  jacobianDeterminant,bcParams,gasParams,rhoBuffer,rhoVBuffer, &
23  rhoEBuffer,viscousFluxBuffer,numscalar,scalarBuffer, &
24  rhoRHS,rhoVRHS,rhoERHS,scalarRHS, &
25  rhoTarget,rhoVTarget,rhoETarget,scalarTarget)
26 
27  IMPLICIT NONE
28 
29  INTEGER(KIND=4) :: numdim,numscalar,patchnormaldir,gridtype
30  INTEGER(KIND=8) :: numpointsbuffer,buffersizes(numdim)
31  INTEGER(KIND=8) :: patchsizes(numdim),numpointspatch
32  INTEGER(KIND=8) :: numpatchpointsop
33  INTEGER(KIND=8) :: patchpointsop(numpatchpointsop)
34  REAL(KIND=8) :: gridmetric(numdim*numdim*numpointsbuffer)
35  REAL(KIND=8) :: jacobiandeterminant(numpointsbuffer)
36  REAL(KIND=8) :: bcparams(3),gasparams(5)
37  REAL(KIND=8) :: rhobuffer(numpointsbuffer)
38  REAL(KIND=8) :: rhovbuffer(numdim*numpointsbuffer)
39  REAL(KIND=8) :: rhoebuffer(numpointsbuffer)
40  REAL(KIND=8) :: viscousfluxbuffer((numdim+2)*numpointsbuffer)
41  REAL(KIND=8) :: scalarbuffer(numscalar*numpointsbuffer)
42  REAL(KIND=8) :: rhorhs(numpointsbuffer)
43  REAL(KIND=8) :: rhovrhs(numdim*numpointsbuffer)
44  REAL(KIND=8) :: rhoerhs(numpointsbuffer)
45  REAL(KIND=8) :: scalarrhs(numscalar*numpointsbuffer)
46  REAL(KIND=8) :: rhotarget(numpointsbuffer)
47  REAL(KIND=8) :: rhovtarget(numdim*numpointsbuffer)
48  REAL(KIND=8) :: rhoetarget(numpointsbuffer)
49  REAL(KIND=8) :: scalartarget(numscalar*numpointsbuffer)
50 
51  ! gasParams(:) = Cref Gamma Cp Re
52  ! bcParams(:) = sigma1 sigma2
53 
54  ! ... Local variables
55  INTEGER(4) :: nd, nauxvars, jj, idim, idir
56  INTEGER(8) :: nc, l0, ipoint, pointindex, metricoffset
57  INTEGER(4) :: normdir, sgn ! , gas_dv_model
58  REAL(8) :: dsgn, sndspdref2, invtempref, tempref, gamref, spcgasconst
59  REAL(8) :: xi_x, xi_y, xi_z, xi_t, bndry_h,sbpboundaryweight
60  REAL(8) :: xi_x_tilde, xi_y_tilde, xi_z_tilde, gridjacobian
61  REAL(8) :: rhofac, penaltyfac, scalarpenaltyfac, scalardiff
62 
63  REAL(8), POINTER :: xix(:)
64  REAL(8) :: metricdata(9)
65  ! REAL(KIND=8), Pointer :: XI_TAU(:,:), cv(:,:), gv(:,:)
66  REAL(8), POINTER :: pnorm_vec(:) ! , rhs(:,:)
67  REAL(8), POINTER :: lambda(:,:), gi1(:), ub(:), tinv(:,:), tmat(:,:)
68  REAL(8), POINTER :: aprime(:,:), gi2(:), matx(:,:), penalty(:) ,correction(:)
69  REAL(8) :: norm_vec(3), gamma, ucon, ubaux,xh
70  REAL(8) :: sigmai1, sigmai2, sigmaaux, sfactor, ibfac_local
71  REAL(8) :: spd_snd, re, reinv, sat_sigmai2_ff, sat_sigmai1_ff
72 
73  LOGICAL :: triggered
74  triggered = .false.
75 
76 
77  ! ... Problem size
78  nd = numdim
79  nc = numpointsbuffer
80  ! N(1:MAX_ND) = 1; Np(1:MAX_ND) = 1;
81  ! Do j = 1, ND
82  ! N(j) = grid_ie(j) - grid_is(j) + 1
83  ! Np(j) = patch_ie(j) - patch_is(j) + 1
84  ! End Do
85  nauxvars = numscalar
86 
87  ! ... Useful quantities about this boundary
88  normdir = abs(patchnormaldir)
89  sgn = normdir / patchnormaldir
90  dsgn = dble(sgn)
91 
92 
93  ! ... Reference quantities
94  ! gasParams = (sndspdref,gamma,Cp,1/Re)
95  ! MJA changed for new nonDimensionalization
96  !sndspdref2 = gasParams(1)*gasParams(1) ! input%sndspdref * input%sndspdref
97  !invtempref = gasParams(3)/sndspdref2 ! input%Cpref / sndspdref2
98  !tempref = 1.0_8 / invtempref
99  !gamref = gasParams(2)
100  !REinv = gasParams(4)
101  gamref = gasparams(2)
102  spcgasconst= gasparams(3)
103  sndspdref2 = gasparams(5)*gasparams(5)
104  tempref = gasparams(4)
105  invtempref = tempref
106  re = gasparams(1)
107  reinv = 0.0_8
108  IF(re .gt. 0) reinv = 1.0_8/re
109 
110 
111  ! ... BC Constants
112  ! bcParams = (sigma1,sigma2,sbpBoundaryStencilWeight)
113  sat_sigmai1_ff = bcparams(1) ! input%SAT_sigmaI1_FF
114  sat_sigmai2_ff = bcparams(2) ! input%SAT_sigmaI2_FF
115  sbpboundaryweight = bcparams(3)
116 
117  ! gas_dv_model = input%gas_dv_model
118 
119 ! WRITE(*,*) 'NORMAL DIR: ',patchNormalDir
120 ! WRITE(*,*) 'SNDSPDREF2 = ',sndspdref2
121 ! WRITE(*,*) 'tempref = ',tempref
122 ! WRITE(*,*) 'gamref = ',gamref
123 ! WRITE(*,*) 'sigma1 = ',SAT_sigmaI1_FF
124 ! WRITE(*,*) 'sigma2 = ',SAT_sigmaI2_FF
125 ! WRITE(*,*) 'REInv = ',REinv
126 ! WRITE(*,*) 'Bweight = ',sbpBoundaryWeight
127 ! WRITE(*,*) 'JAC = ',jacobianDeterminant
128 
129  ! ... memory allocation
130  ALLOCATE(xix(4))
131  ALLOCATE(gi1(nd+2), ub(nd+2), tmat(nd+2,nd+2), tinv(nd+2,nd+2), penalty(nd+2),correction(nd+2))
132  ALLOCATE(gi2(nd+2), aprime(nd+2,nd+2), lambda(nd+2,nd+2), matx(nd+2,nd+2), pnorm_vec(nd))
133 
134 ! WRITE(*,*) 'Farfield ',sgn,normDir
135 
136  xix(:) = 0.0_8
137  IF(gridtype <= unirect) THEN
138  xix(normdir) = gridmetric(normdir)
139  metricoffset = 0
140  gridjacobian = jacobiandeterminant(1)
141  ELSE IF(gridtype == rectilinear) THEN
142  metricoffset = (normdir-1)*numpointsbuffer
143  ELSE
144  metricoffset = (normdir-1)*numdim*numpointsbuffer
145  ENDIF
146 
147 ! WRITE(*,*) 'gridtype,metricoffset = ',gridType,metricOffset
148 
149  ! ... loop over all of the points
150  DO ipoint = 1, numpatchpointsop
151  l0 = patchpointsop(ipoint) + 1 ! assume the points coming in are from C (0-based)
152  ! ... iblank factor (either 1 or zero)
153  ibfac_local = 1.0_8 ! ibfac(l0)
154 
155  ! ... construct the normalized metrics
156  ! XI_X = MT1(l0,t2Map(normDir,1)); XI_Y = 0.0_8; XI_Z = 0.0_8
157  ! if (ND >= 2) XI_Y = MT1(l0,t2Map(normDir,2))
158  ! if (ND == 3) XI_Z = MT1(l0,t2Map(normDir,3))
159  ! XI_T = XI_TAU(l0,normDir)
160  IF(gridtype == curvilinear ) THEN
161  DO idir = 1,numdim
162  xix(idir) = gridmetric(metricoffset + (idir-1)*numpointsbuffer +l0)
163  END DO
164  gridjacobian = jacobiandeterminant(l0)
165  ELSE IF(gridtype == rectilinear) THEN
166  xix(normdir) = gridmetric(metricoffset + l0)
167  gridjacobian = jacobiandeterminant(l0)
168  ENDIF
169 
170  xi_x = xix(1) ! gridMetric(1)
171  xi_y = xix(2)
172  xi_z = xix(3)
173  xi_t = xix(4)
174 
175  IF(numdim == 3) THEN
176 
177  metricdata(1) = xix(1)
178  metricdata(2) = xix(2)
179  metricdata(3) = xix(3)
180 
181  metricdata(4:9) = 0.0_8
182 
183  IF(gridtype < rectilinear) THEN
184  idir = normdir+1
185  IF(idir > 3) idir = idir - 3
186  metricdata(4+(idir-1)) = gridmetric(idir)
187  idir = idir + 1
188  IF(idir > 3) idir = idir - 3
189  metricdata(7+(idir-1)) = gridmetric(idir)
190  ELSE IF(gridtype == rectilinear) THEN
191  idir = normdir + 1
192  IF(idir > 3) idir = idir - 3
193  metricdata(4+(idir-1)) = gridmetric((idir-1)*numpointsbuffer+l0)
194  idir = idir + 1
195  IF(idir > 3) idir = idir - 3
196  metricdata(7+(idir-1)) = gridmetric((idir-1)*numpointsbuffer+l0)
197  ELSE
198  idir = idir+1
199  IF(idir > 3) idir = idir - 3
200  metricoffset = (idir-1)*numdim*numpointsbuffer
201  metricdata(4) = gridmetric(metricoffset+l0)
202  metricdata(5) = gridmetric(metricoffset+numpointsbuffer+l0)
203  metricdata(6) = gridmetric(metricoffset+2*numpointsbuffer+l0)
204  idir = idir+1
205  IF(idir > 3) idir = idir - 3
206  metricoffset = (idir-1)*numdim*numpointsbuffer
207  metricdata(7) = gridmetric(metricoffset+l0)
208  metricdata(8) = gridmetric(metricoffset+numpointsbuffer+l0)
209  metricdata(9) = gridmetric(metricoffset+2*numpointsbuffer+l0)
210  ENDIF
211  END IF
212 
213  ! IF(TRIGGERED.eqv..FALSE.) THEN
214  ! WRITE(*,*) 'XI_X: ',XI_X
215  ! WRITE(*,*) 'XI_Y: ',XI_Y
216  ! WRITE(*,*) 'XI_Z: ',XI_Z
217  ! ENDIF
218 
219  ! ... multiply by the Jacobian
220  xi_x = xi_x * gridjacobian ! JAC(l0)
221  xi_y = xi_y * gridjacobian ! JAC(l0)
222  xi_z = xi_z * gridjacobian ! JAC(l0)
223  xi_t = xi_t * gridjacobian ! JAC(l0)
224 
225  ! ... normalized metrics
226  xh = sqrt((xi_x*xi_x)+(xi_y*xi_y)+(xi_z*xi_z))
227  bndry_h = 1.0_8 / xh
228  xi_x_tilde = xi_x * bndry_h
229  xi_y_tilde = xi_y * bndry_h
230  xi_z_tilde = xi_z * bndry_h
231 
232 
233  ! ... inviscid penalty parameter
234  sigmai1 = -sat_sigmai1_ff * dsgn
235 
236  ! ... pull off the boundary data
237  ! uB(:) = cv(l0,:)
238  ub(1) = rhobuffer(l0)
239  DO idim = 1, numdim
240  ub(1+idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)
241  END DO
242  ub(numdim+2) = rhoebuffer(l0)
243 
244  ! ... target data
245  ! gI1(:) = cvTarget(l0,:)
246  gi1(1) = rhotarget(l0)
247  DO idim = 1, numdim
248  gi1(1+idim) = rhovtarget(l0+(idim-1)*numpointsbuffer)
249  END DO
250  gi1(numdim+2) = rhoetarget(l0)
251 
252  ! ... follow Magnus's algorithm (norm_vec has unit length)
253  norm_vec(1:3) = (/ xi_x_tilde, xi_y_tilde, xi_z_tilde /)
254  pnorm_vec(:) = norm_vec(1:nd) * xh
255 
256  ! ... use target data to make implicit work better?
257  gamma = gamref ! gv(l0,1)
258  ! If (gas_dv_model == DV_MODEL_IDEALGAS .or. &
259  ! gas_dv_model == DV_MODEL_IDEALGAS_MIXTURE) then
260  ! IF(TRIGGERED.eqv..FALSE.) THEN
261  ! WRITE(*,*) 'XI_X_TILDE: ',XI_X_TILDE
262  ! WRITE(*,*) 'XI_Y_TILDE: ',XI_Y_TILDE
263  ! WRITE(*,*) 'XI_Z_TILDE: ',XI_Z_TILDE
264  ! WRITE(*,*) 'XI_X: ',XI_X
265  ! WRITE(*,*) 'XI_Y: ',XI_Y
266  ! WRITE(*,*) 'XI_Z: ',XI_Z
267  ! WRITE(*,*) 'BOUNDARYH: ',bndry_h
268  ! WRITE(*,*) 'GridJacobian: ',gridJacobian
269  ! WRITE(*,*) 'STATE: ',uB(:)
270  ! WRITE(*,*) 'TARGET:',gI1(:)
271 
272  ! ENDIF
273  ! subroutine SAT_Form_Roe_Matrices_Modified(u_in, u_out, gamma, norm, MT1, dir, ND, tmat, tinv, lambda)
274 
275  CALL sat_form_roe_matrices(nd, ub, gi1, gamma, pnorm_vec, tmat, tinv, lambda)
276 ! CALL SAT_Form_Roe_Matrices_test(ND, uB, gI1, gamma, pnorm_vec, metricData,Tmat, Tinv, Lambda)
277 ! CALL SAT_Form_Roe_Matrices_Modified(uB, gI1, gamma, pnorm_vec, metricData, ND, normDir, Tmat, Tinv, Lambda)
278 
279  ! Else If (gas_dv_model == DV_MODEL_THERMALLY_PERFECT) then
280  !
281  ! Call SAT_Form_Roe_Matrices_TP(uB, gI1, ND, gamma, gamref, pnorm_vec, Tmat, Tinv, Lambda)
282  ! End If
283 
284  ! ... save ucon, speed_sound
285  ucon = lambda(1,1)
286  spd_snd = lambda(nd+1,nd+1) - ucon
287 
288  ! ... only pick those Lambda that are incoming
289  If ( sgn == 1 ) Then
290  Do jj = 1, nd+2
291  lambda(jj,jj) = max(lambda(jj,jj),0.0_8)
292  End Do
293  Else
294  Do jj = 1, nd+2
295  lambda(jj,jj) = min(lambda(jj,jj),0.0_8)
296  End Do
297  End If
298 
299  ! ... modify Lambda if subsonic outflow
300  ! ... left boundary
301  if (sgn == 1 .and. &
302  ucon < 0.0_8 .and. ucon + spd_snd > 0.0_8) then
303  lambda(nd+2,nd+1) = lambda(nd+2,nd+1) - ucon - spd_snd
304  ! right boundary
305  else if (sgn == -1 .and. &
306  ucon > 0.0_8 .and. ucon - spd_snd < 0.0_8) then
307  lambda(nd+1,nd+2) = lambda(nd+1,nd+2) - ucon + spd_snd
308  end if
309 
310  ! ... multiply Lambda' into X
311  matx = matmul(lambda,tinv)
312 
313  ! ... compute the characteristic matrix
314  aprime = matmul(tmat,matx)
315 
316  ! ... subtract off the target
317  ub(1:nd+2) = ub(1:nd+2) - gi1(1:nd+2)
318 
319  ! ... compute the characteristic penalty vector
320  gi2(1:nd+2) = matmul(aprime,ub)
321 
322  penaltyfac = sbpboundaryweight*sigmai1
323  correction(:) = penaltyfac*gi2(:)
324  ! ... compute penalty (Bndry_h is included in Lambda already)
325  ! Do jj = 1, ND+2
326  ! rhs(l0,jj) = rhs(l0,jj) + penaltyFac * gI2(jj)
327  ! End Do
328  rhorhs(l0) = rhorhs(l0) + correction(1)
329  DO idim = 1,numdim
330  pointindex = l0 + (idim-1)*numpointsbuffer
331  rhovrhs(pointindex) = rhovrhs(pointindex) + correction(1+idim)
332  END DO
333  rhoerhs(l0) = rhoerhs(l0) + correction(numdim+2)
334 
335 ! IF(TRIGGERED.eqv..FALSE.) THEN
336 ! WRITE(*,*) 'RHS: ',rhoRHS(l0),rhoVRHS(l0),rhoVRHS(l0+numPointsBuffer),rhoERHS(l0)
337 ! WRITE(*,*) 'CORRECTION: ',correction(:)
338 ! TRIGGERED = .TRUE.
339 ! ENDIF
340 
341  ! ... SAT FARFIELD (species treatment)
342  IF(numscalar > 0) THEN
343  If (dsgn * ucon > 0) Then
344  sfactor = ucon * ibfac_local * penaltyfac
345  rhofac = rhobuffer(l0)/rhotarget(l0)
346  scalarpenaltyfac = sfactor*rhofac
347  Do jj = 1, numscalar
348  pointindex = (jj-1)*numpointsbuffer + l0
349  scalardiff = (scalarbuffer(pointindex) - scalartarget(pointindex))
350  scalarrhs(pointindex) = scalarrhs(pointindex) + scalarpenaltyfac * scalardiff
351  ! rhs_auxVars(l0,jj) = rhs_auxVars(l0,jj) + Sfactor * &
352  ! (auxVars(l0,jj) - auxVarsTarget(l0,jj)/cvTarget(l0,1)*cv(l0,1))
353  End do
354  End If
355  END IF
356 
357  ! VISCOUS PART
358  IF(reinv > 0.0) THEN
359  ! ... factor
360 ! sigmaI2 = dsgn * SAT_sigmaI2_FF * SBP_invPdiag_block1(1) / bndry_h !* JAC(l0)
361  sigmai2 = dsgn * sat_sigmai2_ff*sbpboundaryweight/bndry_h !* JAC(l0)
362 
363  ! ... target state
364  do jj = 1, nd+2
365  gi2(jj) = 0.0_8
366  end do
367 
368  ! ... boundary data (already includes 1/Re factor)
369  ub(:) = viscousfluxbuffer(l0) * xi_x_tilde
370  if (nd >= 2) &
371  ub(:) = ub(:) + viscousfluxbuffer(l0) * xi_y_tilde
372  if (nd == 3) &
373  ub(:) = ub(:) + viscousfluxbuffer(l0) * xi_z_tilde
374  ub(:) = 0
375  ! ... penalty term
376  penalty(:) = sigmai2 * (ub(:) - gi2(:))
377 
378  ! ... add to the rhs
379  rhorhs(l0) = rhorhs(l0) + ibfac_local * penalty(1)
380  DO idim = 1,numdim
381  pointindex = l0 + (idim-1)*numpointsbuffer
382  rhovrhs(pointindex) = rhovrhs(pointindex) + ibfac_local*penalty(1+idim)
383  END DO
384  rhoerhs(l0) = rhoerhs(l0) + ibfac_local * penalty(numdim+2)
385 
386  ! species transport
387 ! sigmaAux = ibfac_local * sigmaI2
388 ! do jj = 1, nAuxVars
389 ! ubAux = patch_ViscousFluxAux(jj,lp,normDir)
390 ! rhs_auxVars(l0,jj) = rhs_auxVars(l0,jj) + sigmaAux* ubAux
391 ! end do
392 
393  END IF
394  END DO
395 
396  DEALLOCATE(xix)
397  DEALLOCATE(correction)
398  DEALLOCATE(gi1, ub, tmat, tinv, penalty)
399  DEALLOCATE(gi2, aprime, lambda, matx, pnorm_vec)
400 
401  RETURN
402 
403  END SUBROUTINE farfield
404 
405  SUBROUTINE slip_adiabatic( &
406  numDim,bufferSizes,numPointsBuffer,patchNormalDir,patchSizes, &
407  numPointsPatch,numPatchPointsOp,patchPointsOp,gridType, &
408  gridMetric,jacobianDeterminant,bcParams,gasParams, &
409  rhoBuffer,rhoVBuffer,rhoEBuffer,numscalar,scalarBuffer, &
410  rhoRHS,rhoVRHS,rhoERHS, scalarRHS,rhoTarget,rhoVTarget, &
411  rhoETarget,scalarTarget)
412 
413  IMPLICIT NONE
414 
415  INTEGER(KIND=4) :: numdim,numscalar,patchnormaldir,gridtype
416  INTEGER(KIND=8) :: numpointsbuffer,buffersizes(numdim)
417  INTEGER(KIND=8) :: patchsizes(numdim),numpointspatch
418  INTEGER(KIND=8) :: numpatchpointsop
419  INTEGER(KIND=8) :: patchpointsop(numpatchpointsop)
420  REAL(KIND=8) :: gridmetric(numdim*numdim*numpointsbuffer)
421  REAL(KIND=8) :: jacobiandeterminant(numpointsbuffer)
422  REAL(KIND=8) :: bcparams(3),gasparams(5)
423  REAL(KIND=8) :: rhobuffer(numpointsbuffer)
424  REAL(KIND=8) :: rhovbuffer(numdim*numpointsbuffer)
425  REAL(KIND=8) :: rhoebuffer(numpointsbuffer)
426  REAL(KIND=8) :: scalarbuffer(numscalar*numpointsbuffer)
427  REAL(KIND=8) :: rhorhs(numpointsbuffer)
428  REAL(KIND=8) :: rhovrhs(numdim*numpointsbuffer)
429  REAL(KIND=8) :: rhoerhs(numpointsbuffer)
430  REAL(KIND=8) :: scalarrhs(numscalar*numpointsbuffer)
431  REAL(KIND=8) :: rhotarget(numpointsbuffer)
432  REAL(KIND=8) :: rhovtarget(numdim*numpointsbuffer)
433  REAL(KIND=8) :: rhoetarget(numpointsbuffer)
434  REAL(KIND=8) :: scalartarget(numscalar*numpointsbuffer)
435 
436 
437  ! ... Local variables
438  INTEGER(4) :: nd, nauxvars, jj, idim, idir
439  INTEGER(8) :: nc, l0, ipoint, pointindex,metricoffset
440  INTEGER(4) :: normdir, sgn ! , gas_dv_model
441  ! INTEGER, Dimension(MAX_ND) :: N, Np
442  ! INTEGER, Pointer :: iblank(:), t2Map(:,:)
443  ! INTEGER, Pointer :: grid_is(:), grid_ie(:), patch_is(:), patch_ie(:)
444  ! TYPE (t_grid), Pointer :: grid
445  ! TYPE (t_mixt), Pointer :: state
446  ! TYPE (t_mixt_input), Pointer :: input
447  REAL(8) :: dsgn, sndspdref2, invtempref, tempref, gamref, spcgasconst
448  REAL(8) :: xi_x, xi_y, xi_z, xi_t, bndry_h,sbpboundaryweight
449  REAL(8) :: xi_x_tilde, xi_y_tilde, xi_z_tilde, gridjacobian
450  REAL(8) :: rhofac, penaltyfac, scalarpenaltyfac, scalardiff
451 
452  REAL(8), POINTER :: xix(:), bvelocity(:), vwall(:), vwalltarget(:), vtarget(:)
453  REAL(8) :: metricdata(9)
454  ! REAL(KIND=8), Pointer :: XI_TAU(:,:), cv(:,:), gv(:,:)
455  REAL(8), POINTER :: pnorm_vec(:) ! , rhs(:,:)
456  ! REAL(KIND=8), Pointer :: patch_ViscousFlux(:,:,:)
457  ! REAL(KIND=8), Pointer :: patch_ViscousFluxAux(:,:,:)
458  ! REAL(KIND=8), Pointer :: patch_ViscousRHSAux(:,:,:), rhs_auxVars(:,:)
459  ! REAL(KIND=8), Pointer :: auxVars(:,:), Lambda(:,:), ibfac(:)
460  REAL(8), POINTER :: lambda(:,:), gi1(:), ub(:), tinv(:,:), tmat(:,:)
461  ! REAL(KIND=8), Pointer :: gI1(:), ub(:), Tinv(:,:), Tmat(:,:)
462  REAL(8), POINTER :: aprime(:,:), gi2(:), matx(:,:), penalty(:) ,correction(:)
463  ! REAL(KIND=8), Pointer :: auxVarsTarget(:,:), cvTarget(:,:)
464  ! REAL(KIND=8), Pointer :: SBP_invPdiag_block1(:), JAC(:), MT1(:,:)
465  REAL(8) :: norm_vec(3), gamma, ucon, ubaux, specrad, vdotxi, xh
466  REAL(8) :: sigmai1, sigmai2, sigmaaux, sfactor, ibfac_local
467  REAL(8) :: spd_snd, re, reinv, sat_sigmai1,pressuretarget
468 
469  LOGICAL :: triggered
470  triggered = .false.
471 
472  ! ... Main data containers
473  ! patch_gridID = patch%gridID
474  ! grid => region%grid(patch_gridID)
475  ! state => region%state(patch_gridID)
476  ! input => grid%input
477 
478  ! ... Pointer dereferences
479  ! grid_is => grid%is
480  ! grid_ie => grid%ie
481  ! patch_is => patch%is
482  ! patch_ie => patch%ie
483  ! cv => state%cv
484  ! gv => state%gv
485  ! cvTarget => state%cvTarget
486  ! rhs => state%rhs
487  ! iblank => grid%iblank
488  ! JAC => grid%JAC
489  ! MT1 => grid%MT1
490  ! ibfac => grid%ibfac
491  ! t2Map => region%global%t2Map
492  ! XI_TAU => grid%XI_TAU
493  ! rhs_auxVars => state%rhs_auxVars
494  ! patch_ViscousFlux => patch%ViscousFlux
495  ! patch_ViscousFluxAux => patch%ViscousFluxAux
496  ! SBP_invPdiag_block1 => grid%SBP_invPdiag_block1
497  ! auxVars => state%auxVars
498  ! auxVarsTarget => state%auxVarsTarget
499 
500  ! ... Problem size
501  nd = numdim
502  nc = numpointsbuffer
503 
504  nauxvars = numscalar
505 
506  ! ... Useful quantities about this boundary
507  normdir = abs(patchnormaldir)
508  sgn = normdir / patchnormaldir
509  dsgn = dble(sgn)
510 
511 
512  ! ... Reference quantities
513  ! gasParams = (sndspdref,gamma,Cp,1/Re)
514  ! bcParams(:) = sigma1 sigma2
515  ! MJA changed for new nonDimensionalization
516  !sndspdref2 = gasParams(1)*gasParams(1) ! input%sndspdref * input%sndspdref
517  !invtempref = gasParams(3)/sndspdref2 ! input%Cpref / sndspdref2
518  !tempref = 1.0_8 / invtempref
519  !gamref = gasParams(2)
520  !REinv = gasParams(4)
521  gamref = gasparams(2)
522  spcgasconst= gasparams(3)
523  sndspdref2 = gasparams(5)*gasparams(5)
524  tempref = gasparams(4)
525  invtempref = tempref
526  re = gasparams(1)
527  reinv = 1.0_8/re
528 
529  ! ... BC Constants
530  ! bcParams = (sigma1,sbpBoundaryStencilWeight)
531  sat_sigmai1 = bcparams(1) ! input%SAT_sigmaI1_FF
532  sbpboundaryweight = bcparams(2)
533 
534  ! gas_dv_model = input%gas_dv_model
535 
536  ! ... memory allocation
537  ALLOCATE(xix(4),vwall(nd),bvelocity(nd),vwalltarget(nd),vtarget(nd))
538  ALLOCATE(gi1(nd+2), ub(nd+2), tmat(nd+2,nd+2), tinv(nd+2,nd+2), penalty(nd+2),correction(nd+2))
539  ALLOCATE(gi2(nd+2), aprime(nd+2,nd+2), lambda(nd+2,nd+2), matx(nd+2,nd+2), pnorm_vec(nd))
540 
541  xix(:) = 0.0_8
542  IF(gridtype <= unirect) THEN
543  xix(normdir) = gridmetric(normdir)
544  metricoffset = 0
545  gridjacobian = jacobiandeterminant(1)
546  ELSE IF(gridtype == rectilinear) THEN
547  metricoffset = (normdir-1)*numpointsbuffer
548  ELSE
549  metricoffset = (normdir-1)*numdim*numpointsbuffer
550  ENDIF
551 
552 ! XIX(:) = 0.0_8
553 ! XIX(normDir) = gridMetric(normDir) ! only diagonal terms are non-zero
554 
555 ! WRITE(*,*) 'XIX = ',XIX
556 ! WRITE(*,*) 'NORMAL DIR: ',patchNormalDir
557 ! WRITE(*,*) 'SNDSPDREF2 = ',sndspdref2
558 ! WRITE(*,*) 'tempref = ',tempref
559 ! WRITE(*,*) 'gamref = ',gamref
560 ! WRITE(*,*) 'sigma1 = ',SAT_sigmaI1
561 ! WRITE(*,*) 'REInv = ',REinv
562 ! WRITE(*,*) 'Bweight = ',sbpBoundaryWeight
563 ! WRITE(*,*) 'JAC = ',jacobianDeterminant
564 
565 ! TRIGGERED = .TRUE.
566 
567  ! ... loop over all of the points
568  DO ipoint = 1, numpatchpointsop
569 
570  l0 = patchpointsop(ipoint) + 1 ! assume the points coming in are from C (0-based)
571  ! lp = (k-patch_is(3))*Np(1)*Np(2) + (j-patch_is(2))*Np(1) + (i-patch_is(1)+1)
572 
573  ! IF(iPoint == (numPatchPointsOp/2)) TRIGGERED = .FALSE.
574 
575  ! ... iblank factor (either 1 or zero)
576  ibfac_local = 1.0_8 ! ibfac(l0)
577 
578  ! ... construct the normalized metrics
579  ! XI_X = MT1(l0,t2Map(normDir,1)); XI_Y = 0.0_8; XI_Z = 0.0_8
580  ! if (ND >= 2) XI_Y = MT1(l0,t2Map(normDir,2))
581  ! if (ND == 3) XI_Z = MT1(l0,t2Map(normDir,3))
582  ! XI_T = XI_TAU(l0,normDir)
583  ! IF(ND > 1) XI_Y = gridMetric(2)
584  ! IF(ND > 2) XI_Z = gridMetric(3)
585 
586 
587  IF(gridtype == curvilinear ) THEN
588  DO idir = 1,numdim
589  xix(idir) = gridmetric(metricoffset + (idir-1)*numpointsbuffer +l0)
590  END DO
591  gridjacobian = jacobiandeterminant(l0)
592  ELSE IF(gridtype == rectilinear) THEN
593  xix(normdir) = gridmetric(metricoffset + l0)
594  gridjacobian = jacobiandeterminant(l0)
595  ENDIF
596 
597  ! ... multiply by the Jacobian
598  xi_x = xix(1) * gridjacobian ! JAC(l0)
599  xi_y = xix(2) * gridjacobian ! JAC(l0)
600  xi_z = xix(3) * gridjacobian ! JAC(l0)
601  xi_t = xix(4) * gridjacobian ! JAC(l0)
602 
603  IF(numdim == 3) THEN
604 
605  metricdata(1) = xix(1)
606  metricdata(2) = xix(2)
607  metricdata(3) = xix(3)
608  metricdata(4:9) = 0.0_8
609 
610  IF(gridtype < rectilinear) THEN
611  idir = normdir+1
612  IF(idir > 3) idir = idir - 3
613  metricdata(4+(idir-1)) = gridmetric(idir)
614  idir = idir + 1
615  IF(idir > 3) idir = idir - 3
616  metricdata(7+(idir-1)) = gridmetric(idir)
617  ELSE IF(gridtype == rectilinear) THEN
618  idir = normdir + 1
619  IF(idir > 3) idir = idir - 3
620  metricdata(4+(idir-1)) = gridmetric((idir-1)*numpointsbuffer+l0)
621  idir = idir + 1
622  IF(idir > 3) idir = idir - 3
623  metricdata(7+(idir-1)) = gridmetric((idir-1)*numpointsbuffer+l0)
624  ELSE
625  idir = idir+1
626  IF(idir > 3) idir = idir - 3
627  metricoffset = (idir-1)*numdim*numpointsbuffer
628  metricdata(4) = gridmetric(metricoffset+l0)
629  metricdata(5) = gridmetric(metricoffset+numpointsbuffer+l0)
630  metricdata(6) = gridmetric(metricoffset+2*numpointsbuffer+l0)
631  idir = idir+1
632  IF(idir > 3) idir = idir - 3
633  metricoffset = (idir-1)*numdim*numpointsbuffer
634  metricdata(7) = gridmetric(metricoffset+l0)
635  metricdata(8) = gridmetric(metricoffset+numpointsbuffer+l0)
636  metricdata(9) = gridmetric(metricoffset+2*numpointsbuffer+l0)
637  ENDIF
638  END IF
639 
640  ! ... normalized metrics
641  xh = sqrt( (xi_x * xi_x) + (xi_y * xi_y) + (xi_z * xi_z) )
642  bndry_h = 1.0_8 / xh
643  xi_x_tilde = xi_x * bndry_h
644  xi_y_tilde = xi_y * bndry_h
645  xi_z_tilde = xi_z * bndry_h
646 
647 
648  ! ... follow Magnus's algorithm (norm_vec has unit length)
649  norm_vec(1:3) = (/ xi_x_tilde, xi_y_tilde, xi_z_tilde /)
650  pnorm_vec(:) = norm_vec(1:nd) * xh
651 
652 
653  ! ... penalty parameter
654  sigmai1 = -sat_sigmai1 * dsgn
655 
656  ! ... inviscid penalty parameter
657  ! sigmaI1 = -SAT_sigmaI1_FF * dsgn
658 
659  ! ... pull off the boundary data
660  ! uB(:) = cv(l0,:)
661  ub(1) = rhobuffer(l0)
662  vdotxi = 0
663  DO idim = 1, numdim
664  ub(1+idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)
665  bvelocity(idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)/rhobuffer(l0)
666  vdotxi = vdotxi + bvelocity(idim)*xix(idim)*gridjacobian*bndry_h
667  END DO
668  ub(numdim+2) = rhoebuffer(l0)
669 
670  vwall(1:numdim) = vdotxi*norm_vec(1:numdim)
671  ! NON MOVING GRIDS!!
672  vwalltarget(1:numdim) = 0.0_8
673  ! = dot_product(XYZ_TAU(l0,1:ND),norm_vec(1:ND)) * norm_vec(1:ND)
674 
675 
676  DO idim = 1,numdim
677  vtarget(idim) = bvelocity(idim) - ( vwall(idim) - vwalltarget(idim) )
678  END DO
679 ! pressureTarget = pressureBuffer(l0)
680 
681  ! ... target data
682  ! gI1(:) = cvTarget(l0,:)
683  gi1(1) = ub(1)
684  gi1(numdim+2) = ub(numdim+2)
685  DO idim = 1, numdim
686  gi1(1+idim) = gi1(1) * vtarget(idim)
687  gi1(numdim+2) = gi1(numdim+2) + 0.5_8*gi1(1)*(vtarget(idim)**2 - bvelocity(idim)**2)
688  END DO
689 
690 
691  ! ... use target data to make implicit work better?
692  gamma = gamref ! gv(l0,1)
693  ! If (gas_dv_model == DV_MODEL_IDEALGAS .or. &
694  ! gas_dv_model == DV_MODEL_IDEALGAS_MIXTURE) then
695 
696 
697 ! IF(TRIGGERED .EQV. .FALSE.) THEN
698 ! WRITE(*,*) 'L0 = ',l0
699 ! WRITE(*,*) 'XIX = ',XIX
700 ! WRITE(*,*) '1/xh = ',1.0_8/xh
701 ! WRITE(*,*) 'VdotXi = ',vdotxi
702 ! WRITE(*,*) 'bVelocity = ',bVelocity
703 ! WRITE(*,*) 'pnorm: ',pnorm_vec
704 ! WRITE(*,*) 'State: ',uB
705 ! WRITE(*,*) 'Target: ',gI1
706 ! TRIGGERED = .TRUE.
707 ! ENDIF
708 
709  CALL sat_form_roe_matrices(nd, ub, gi1, gamma, pnorm_vec, tmat, tinv, lambda)
710 
711  specrad = 0.0_8
712  ucon = vdotxi*xh
713 
714  ! ... only pick those Lambda that are incoming
715  If ( sgn == 1 ) Then
716  Do jj = 1, nd+2
717  lambda(jj,jj) = max(lambda(jj,jj),0.0_8)
718  specrad = max(abs(lambda(jj,jj)),abs(specrad))
719  End Do
720  Else
721  Do jj = 1, nd+2
722  lambda(jj,jj) = min(lambda(jj,jj),0.0_8)
723  specrad = max(abs(lambda(jj,jj)),abs(specrad))
724  End Do
725  End If
726 
727 ! IF(TRIGGERED.eqv..FALSE.) THEN
728 ! WRITE(*,*) 'XI_X_TILDE: ',XI_X_TILDE
729 ! WRITE(*,*) 'XI_Y_TILDE: ',XI_Y_TILDE
730 ! WRITE(*,*) 'XI_Z_TILDE: ',XI_Z_TILDE
731 ! WRITE(*,*) 'XI_X: ',XI_X
732 ! WRITE(*,*) 'XI_Y: ',XI_Y
733 ! WRITE(*,*) 'XI_Z: ',XI_Z
734 ! WRITE(*,*) 'BOUNDARYH: ',bndry_h
735 ! WRITE(*,*) 'STATE: ',uB(:)
736 ! WRITE(*,*) 'TARGET:',gI1(:)
737 ! WRITE(*,*) 'UCON: ',uCon
738 ! WRITE(*,*) 'SpecRad: ',specRad
739 ! WRITE(*,*) 'vDotXi: ',vDotXi
740 ! ENDIF
741 
742  matx = matmul(lambda,tinv)
743 
744  ! ... compute the characteristic matrix
745  aprime = matmul(tmat,matx)
746 
747  ! ... subtract off the target
748  ub(1:nd+2) = ub(1:nd+2) - gi1(1:nd+2)
749 
750  ! ... compute the characteristic penalty vector
751  gi2(1:nd+2) = matmul(aprime,ub)
752 
753  penaltyfac = sbpboundaryweight*sigmai1
754  correction(:) = penaltyfac*gi2(:)
755 
756  ! IBFAC would be considered here
757  rhorhs(l0) = rhorhs(l0) + ibfac_local * correction(1)
758  DO idim = 1,numdim
759  pointindex = l0 + (idim-1)*numpointsbuffer
760  rhovrhs(pointindex) = rhovrhs(pointindex) + ibfac_local*correction(1+idim)
761  END DO
762  rhoerhs(l0) = rhoerhs(l0) + ibfac_local * correction(numdim+2)
763 
764 ! IF(TRIGGERED.eqv..FALSE.) THEN
765 ! WRITE(*,*) 'CORRECTION(',l0,'): ',correction(:)
766 ! TRIGGERED = .TRUE.
767 ! ENDIF
768 
769  ! ... SAT_SLIP_ADIABATIC (species treatment)
770  IF(numscalar > 0) THEN
771  If (dsgn * ucon > 0) Then
772  sfactor = -ucon * dsgn * ibfac_local * sbpboundaryweight
773  Do jj = 1, numscalar
774  pointindex = (jj-1)*numpointsbuffer + l0
775  scalarrhs(pointindex) = scalarrhs(pointindex) + sfactor*scalarbuffer(pointindex)
776  End do
777  End If
778  END IF
779 
781  ! VISCOUS PART OMITTED (for now)
782 
783 
784 ! if ( REinv > 0.0_8 ) then
785 ! sigmaI2 = ibfac_local * dsgn * 1_8 * SBP_invPdiag_block1(1) * JAC(l0)
786 ! gI2 = 0_8
787 ! uB(:) = patch_ViscousFlux(:,lp,normDir)
788 ! penalty(:) = sigmaI2 * (uB(:) - gI2(:))
789 ! penalty(normDir+1) = 0_8 !normal momentum
790 ! do jj = 1, ND+2
791 ! rhs(l0,jj) = rhs(l0,jj) + penalty(jj)
792 ! end do
793 ! do jj = 1, nAuxVars
794 ! ubAux = patch_ViscousFluxAux(jj,lp,normDir)
795 ! rhs_auxVars(l0,jj) = rhs_auxVars(l0,jj) + sigmaI2* ubAux
796 ! end do
797 ! endif !REinv > 0.0_8
798 
799 
800  END DO
801 
802  DEALLOCATE(xix,vwall,bvelocity,vwalltarget,vtarget)
803  DEALLOCATE(gi1, ub, tmat, tinv, penalty, correction)
804  DEALLOCATE(gi2, aprime, lambda, matx, pnorm_vec)
805 
806  RETURN
807 
808  END SUBROUTINE slip_adiabatic
809 
810  SUBROUTINE noslip_isothermal( &
811  numDim,bufferSizes,numPointsBuffer,patchNormalDir,patchSizes, &
812  numPointsPatch,numPatchPointsOp,patchPointsOp,gridType, &
813  gridMetric,jacobianDeterminant,bcParams,gasParams, &
814  rhoBuffer,rhoVBuffer,rhoEBuffer,numscalar,scalarBuffer, &
815  rhoRHS,rhoVRHS,rhoERHS, scalarRHS,rhoTarget,rhoVTarget, &
816  rhoETarget,scalarTarget,muBuffer,lambdaBuffer)
817 
818  IMPLICIT NONE
819 
820  INTEGER(KIND=4) :: numdim,numscalar,patchnormaldir,gridtype
821  INTEGER(KIND=8) :: numpointsbuffer,buffersizes(numdim)
822  INTEGER(KIND=8) :: patchsizes(numdim),numpointspatch
823  INTEGER(KIND=8) :: numpatchpointsop
824  INTEGER(KIND=8) :: patchpointsop(numpatchpointsop)
825  REAL(KIND=8) :: gridmetric(numdim*numdim*numpointsbuffer)
826  REAL(KIND=8) :: jacobiandeterminant(numpointsbuffer)
827  REAL(KIND=8) :: bcparams(4),gasparams(6)
828  REAL(KIND=8) :: rhobuffer(numpointsbuffer)
829  REAL(KIND=8) :: rhovbuffer(numdim*numpointsbuffer)
830  REAL(KIND=8) :: rhoebuffer(numpointsbuffer)
831  REAL(KIND=8) :: scalarbuffer(numscalar*numpointsbuffer)
832  REAL(KIND=8) :: rhorhs(numpointsbuffer)
833  REAL(KIND=8) :: rhovrhs(numdim*numpointsbuffer)
834  REAL(KIND=8) :: rhoerhs(numpointsbuffer)
835  REAL(KIND=8) :: scalarrhs(numscalar*numpointsbuffer)
836  REAL(KIND=8) :: rhotarget(numpointsbuffer)
837  REAL(KIND=8) :: rhovtarget(numdim*numpointsbuffer)
838  REAL(KIND=8) :: rhoetarget(numpointsbuffer)
839  REAL(KIND=8) :: scalartarget(numscalar*numpointsbuffer)
840  REAL(KIND=8) :: mubuffer(numpointsbuffer)
841  REAL(KIND=8) :: lambdabuffer(numpointsbuffer)
842 
843 
844  ! ... Local variables
845  INTEGER(4) :: nd, nauxvars, jj, idim, idir
846  INTEGER(8) :: nc, l0, ipoint, pointindex, metricoffset
847  INTEGER(4) :: normdir, sgn ! , gas_dv_model
848  INTEGER(4) :: numdebugpoints, idebugpoint,debugpoints(6)
849  ! INTEGER, Dimension(MAX_ND) :: N, Np
850  ! INTEGER, Pointer :: iblank(:), t2Map(:,:)
851  ! INTEGER, Pointer :: grid_is(:), grid_ie(:), patch_is(:), patch_ie(:)
852  ! TYPE (t_grid), Pointer :: grid
853  ! TYPE (t_mixt), Pointer :: state
854  ! TYPE (t_mixt_input), Pointer :: input
855  REAL(8) :: dsgn, sndspdref2, invtempref, tempref, gamref, spcgasconst
856  REAL(8) :: xi_x, xi_y, xi_z, xi_t, bndry_h,sbpboundaryweight, bcwalltemp
857  REAL(8) :: xi_x_tilde, xi_y_tilde, xi_z_tilde, gridjacobian
858  REAL(8) :: rhofac, penaltyfac, scalarpenaltyfac, scalardiff
859 
860  REAL(8), POINTER :: xix(:), bvelocity(:), vwall(:), vwalltarget(:), vtarget(:)
861  REAL(8) :: metricdata(9)
862  ! REAL(KIND=8), Pointer :: XI_TAU(:,:), cv(:,:), gv(:,:)
863  REAL(8), POINTER :: pnorm_vec(:) ! , rhs(:,:)
864  ! REAL(KIND=8), Pointer :: patch_ViscousFlux(:,:,:)
865  ! REAL(KIND=8), Pointer :: patch_ViscousFluxAux(:,:,:)
866  ! REAL(KIND=8), Pointer :: patch_ViscousRHSAux(:,:,:), rhs_auxVars(:,:)
867  ! REAL(KIND=8), Pointer :: auxVars(:,:), Lambda(:,:), ibfac(:)
868  REAL(8), POINTER :: lambda(:,:), gi1(:), ub(:), tinv(:,:), tmat(:,:)
869  ! REAL(KIND=8), Pointer :: gI1(:), ub(:), Tinv(:,:), Tmat(:,:)
870  REAL(8), POINTER :: aprime(:,:), gi2(:), matx(:,:), penalty(:) ,correction(:)
871  ! REAL(KIND=8), Pointer :: auxVarsTarget(:,:), cvTarget(:,:)
872  ! REAL(KIND=8), Pointer :: SBP_invPdiag_block1(:), JAC(:), MT1(:,:)
873  REAL(8) :: norm_vec(3), gamma, ucon, ubaux, specrad, vdotxi, xh
874  REAL(8) :: sigmai1, sigmai2, sigmaaux, sfactor, ibfac_local, t_wall
875  REAL(8) :: spd_snd, re, reinv, sat_sigmai1, sat_sigmai2, pressuretarget
876  REAL(8) :: viscouspenaltyscale
877 
878  LOGICAL :: triggered
879  triggered = .false.
880 
881  ! ... Main data containers
882  ! patch_gridID = patch%gridID
883  ! grid => region%grid(patch_gridID)
884  ! state => region%state(patch_gridID)
885  ! input => grid%input
886 
887  ! ... Pointer dereferences
888  ! grid_is => grid%is
889  ! grid_ie => grid%ie
890  ! patch_is => patch%is
891  ! patch_ie => patch%ie
892  ! cv => state%cv
893  ! gv => state%gv
894  ! cvTarget => state%cvTarget
895  ! rhs => state%rhs
896  ! iblank => grid%iblank
897  ! JAC => grid%JAC
898  ! MT1 => grid%MT1
899  ! ibfac => grid%ibfac
900  ! t2Map => region%global%t2Map
901  ! XI_TAU => grid%XI_TAU
902  ! rhs_auxVars => state%rhs_auxVars
903  ! patch_ViscousFlux => patch%ViscousFlux
904  ! patch_ViscousFluxAux => patch%ViscousFluxAux
905  ! SBP_invPdiag_block1 => grid%SBP_invPdiag_block1
906  ! auxVars => state%auxVars
907  ! auxVarsTarget => state%auxVarsTarget
908 
909  ! ... Problem size
910  nd = numdim
911  nc = numpointsbuffer
912 
913  nauxvars = numscalar
914 
915  ! ... Useful quantities about this boundary
916  normdir = abs(patchnormaldir)
917  sgn = normdir / patchnormaldir
918  dsgn = dble(sgn)
919  numdebugpoints = 0
920 
921  ! ... Reference quantities
922  ! gasParams = (sndspdref,gamma,Cp,1/Re)
923  ! bcParams(:) = sigma1 sigma2
924  ! MJA changed for new nonDimensionalization
925  !sndspdref2 = gasParams(1)*gasParams(1) ! input%sndspdref * input%sndspdref
926  !invtempref = gasParams(3)/sndspdref2 ! input%Cpref / sndspdref2
927  !tempref = 1.0_8 / invtempref
928  !gamref = gasParams(2)
929  !REinv = gasParams(4)
930  gamref = gasparams(2)
931  spcgasconst= gasparams(3)
932  sndspdref2 = gasparams(5)*gasparams(5)
933  tempref = gasparams(4)
934  invtempref = tempref
935  re = gasparams(1)
936  reinv = 1.0_8/re
937  viscouspenaltyscale = gasparams(6)
938 
939  ! ... BC Constants
940  ! bcParams = (sigma1,sbpBoundaryStencilWeight)
941  sat_sigmai1 = bcparams(1) ! input%SAT_sigmaI1_FF
942  sat_sigmai2 = bcparams(2) ! /input%SAT_sigmaI2
943  bcwalltemp = bcparams(3) ! input%bcic_wallTemp
944  sbpboundaryweight = bcparams(4) !
945 
946  !write(*,*) "SAT_sigmaI1 SAT_sigmaI2 bcWallTemp sbpBoundaryWeight tempRef"
947  !write(*,*) SAT_sigmaI1, SAT_sigmaI2, bcWallTemp, sbpBoundaryWeight, tempRef
948 
949  ! gas_dv_model = input%gas_dv_model
950 
951  ! ... memory allocation
952  ALLOCATE(xix(4),vwall(nd),bvelocity(nd),vwalltarget(nd),vtarget(nd))
953  ALLOCATE(gi1(nd+2), ub(nd+2), tmat(nd+2,nd+2), tinv(nd+2,nd+2), penalty(nd+2),correction(nd+2))
954  ALLOCATE(gi2(nd+2), aprime(nd+2,nd+2), lambda(nd+2,nd+2), matx(nd+2,nd+2), pnorm_vec(nd))
955 
956  xix(:) = 0.0_8
957  IF(gridtype <= unirect) THEN
958  xix(normdir) = gridmetric(normdir)
959  metricoffset = 0
960  gridjacobian = jacobiandeterminant(1)
961  ELSE IF(gridtype == rectilinear) THEN
962  metricoffset = (normdir-1)*numpointsbuffer
963  ELSE
964  metricoffset = (normdir-1)*numdim*numpointsbuffer
965  ENDIF
966 
967 ! WRITE(*,*) 'XIX = ',XIX
968 ! WRITE(*,*) 'NORMAL DIR: ',patchNormalDir
969 ! WRITE(*,*) 'SNDSPDREF2 = ',sndspdref2
970 ! WRITE(*,*) 'tempref = ',tempref
971 ! WRITE(*,*) 'gamref = ',gamref
972 ! WRITE(*,*) 'sigma1 = ',SAT_sigmaI1
973 ! WRITE(*,*) 'REInv = ',REinv
974 ! WRITE(*,*) 'Bweight = ',sbpBoundaryWeight
975 ! WRITE(*,*) 'JAC = ',jacobianDeterminant
976 
977 ! TRIGGERED = .TRUE.
978  debugpoints(:) = 0
979  numdebugpoints = 6
980  debugpoints(1) = 91427
981  debugpoints(2) = 91428
982  debugpoints(3) = 91678
983  debugpoints(4) = 114268
984  debugpoints(5) = 114269
985  debugpoints(6) = 114018
986 
987  ! ... loop over all of the points
988  DO ipoint = 1, numpatchpointsop
989 
990  l0 = patchpointsop(ipoint) + 1 ! assume the points coming in are from C (0-based)
991  ! lp = (k-patch_is(3))*Np(1)*Np(2) + (j-patch_is(2))*Np(1) + (i-patch_is(1)+1)
992 
993  ! IF(iPoint == (numPatchPointsOp/2)) TRIGGERED = .FALSE.
994 
995  ! ... iblank factor (either 1 or zero)
996  ibfac_local = 1.0_8 ! ibfac(l0)
997 
998 
999  ! ... construct the normalized metrics
1000  ! XI_X = MT1(l0,t2Map(normDir,1)); XI_Y = 0.0_8; XI_Z = 0.0_8
1001  ! if (ND >= 2) XI_Y = MT1(l0,t2Map(normDir,2))
1002  ! if (ND == 3) XI_Z = MT1(l0,t2Map(normDir,3))
1003  ! XI_T = XI_TAU(l0,normDir)
1004  ! IF(ND > 1) XI_Y = gridMetric(2)
1005  ! IF(ND > 2) XI_Z = gridMetric(3)
1006  IF(gridtype == curvilinear ) THEN
1007  DO idir = 1,numdim
1008  xix(idir) = gridmetric(metricoffset + (idir-1)*numpointsbuffer +l0)
1009  END DO
1010  gridjacobian = jacobiandeterminant(l0)
1011  ELSE IF(gridtype == rectilinear) THEN
1012  xix(normdir) = gridmetric(metricoffset + l0)
1013  gridjacobian = jacobiandeterminant(l0)
1014  ENDIF
1015 
1016 
1017  IF(numdim == 3) THEN
1018 
1019  metricdata(1) = xix(1)
1020  metricdata(2) = xix(2)
1021  metricdata(3) = xix(3)
1022  metricdata(4:9) = 0.0_8
1023 
1024  IF(gridtype < rectilinear) THEN
1025  idir = normdir+1
1026  IF(idir > 3) idir = idir - 3
1027  metricdata(4+(idir-1)) = gridmetric(idir)
1028  idir = idir + 1
1029  IF(idir > 3) idir = idir - 3
1030  metricdata(7+(idir-1)) = gridmetric(idir)
1031  ELSE IF(gridtype == rectilinear) THEN
1032  idir = normdir + 1
1033  IF(idir > 3) idir = idir - 3
1034  metricdata(4+(idir-1)) = gridmetric((idir-1)*numpointsbuffer+l0)
1035  idir = idir + 1
1036  IF(idir > 3) idir = idir - 3
1037  metricdata(7+(idir-1)) = gridmetric((idir-1)*numpointsbuffer+l0)
1038  ELSE
1039  idir = idir+1
1040  IF(idir > 3) idir = idir - 3
1041  metricoffset = (idir-1)*numdim*numpointsbuffer
1042  metricdata(4) = gridmetric(metricoffset+l0)
1043  metricdata(5) = gridmetric(metricoffset+numpointsbuffer+l0)
1044  metricdata(6) = gridmetric(metricoffset+2*numpointsbuffer+l0)
1045  idir = idir+1
1046  IF(idir > 3) idir = idir - 3
1047  metricoffset = (idir-1)*numdim*numpointsbuffer
1048  metricdata(7) = gridmetric(metricoffset+l0)
1049  metricdata(8) = gridmetric(metricoffset+numpointsbuffer+l0)
1050  metricdata(9) = gridmetric(metricoffset+2*numpointsbuffer+l0)
1051  ENDIF
1052  END IF
1053  ! ... multiply by the Jacobian
1054  xi_x = xix(1) * gridjacobian ! JAC(l0)
1055  xi_y = xix(2) * gridjacobian ! JAC(l0)
1056  xi_z = xix(3) * gridjacobian ! JAC(l0)
1057  xi_t = xix(4) * gridjacobian ! JAC(l0)
1058 
1059  ! ... normalized metrics
1060  xh = sqrt( (xi_x * xi_x) + (xi_y * xi_y) + (xi_z * xi_z) )
1061  bndry_h = 1.0_8 / xh
1062  xi_x_tilde = xi_x * bndry_h
1063  xi_y_tilde = xi_y * bndry_h
1064  xi_z_tilde = xi_z * bndry_h
1065 
1066 
1067  ! ... follow Magnus's algorithm (norm_vec has unit length)
1068  norm_vec(1:3) = (/ xi_x_tilde, xi_y_tilde, xi_z_tilde /)
1069  pnorm_vec(:) = norm_vec(1:nd) * xh
1070 
1071  ! ... penalty parameter
1072  sigmai1 = -sat_sigmai1 * dsgn
1073 
1074  ! ... inviscid penalty parameter
1075  ! sigmaI1 = -SAT_sigmaI1_FF * dsgn
1076 
1077  ! ... pull off the boundary data
1078  ! uB(:) = cv(l0,:)
1079  ub(1) = rhobuffer(l0)
1080  vdotxi = 0
1081  DO idim = 1, numdim
1082  ub(1+idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)
1083  bvelocity(idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)/rhobuffer(l0)
1084 ! vDotXi = vDotXi + bVelocity(iDim)*XIX(iDim)*gridJacobian*bndry_h
1085  vdotxi = vdotxi + bvelocity(idim)*norm_vec(idim)
1086  END DO
1087  ub(numdim+2) = rhoebuffer(l0)
1088 
1089  vwall(1:numdim) = vdotxi*norm_vec(1:numdim)
1090  ! NON MOVING GRIDS!!
1091  vwalltarget(1:numdim) = 0.0_8
1092  ! = dot_product(XYZ_TAU(l0,1:ND),norm_vec(1:ND)) * norm_vec(1:ND)
1093 
1094 
1095  DO idim = 1,numdim
1096  vtarget(idim) = bvelocity(idim) - ( vwall(idim) - vwalltarget(idim) )
1097  END DO
1098 
1099 ! pressureTarget = pressureBuffer(l0)
1100 
1101  ! ... target data
1102  ! gI1(:) = cvTarget(l0,:)
1103  gi1(1) = ub(1)
1104  gi1(numdim+2) = ub(numdim+2)
1105  DO idim = 1, numdim
1106  gi1(1+idim) = gi1(1) * vtarget(idim)
1107  gi1(numdim+2) = gi1(numdim+2) + 0.5_8*gi1(1)*(vtarget(idim)**2 - bvelocity(idim)**2)
1108  END DO
1109 
1110 
1111  ! ... use target data to make implicit work better?
1112  gamma = gamref ! gv(l0,1)
1113  ! If (gas_dv_model == DV_MODEL_IDEALGAS .or. &
1114  ! gas_dv_model == DV_MODEL_IDEALGAS_MIXTURE) then
1115 
1116 
1117 ! IF(TRIGGERED .EQV. .FALSE.) THEN
1118 ! WRITE(*,*) 'L0 = ',l0
1119 ! WRITE(*,*) 'XIX = ',XIX
1120 ! WRITE(*,*) '1/xh = ',1.0_8/xh
1121 ! WRITE(*,*) 'VdotXi = ',vdotxi
1122 ! WRITE(*,*) 'bVelocity = ',bVelocity
1123 ! WRITE(*,*) 'pnorm: ',pnorm_vec
1124 ! WRITE(*,*) 'State: ',uB
1125 ! WRITE(*,*) 'Target: ',gI1
1126 ! TRIGGERED = .TRUE.
1127 ! ENDIF
1128 
1129 ! DO iDebugPoint = 1, numDebugPoints
1130 ! IF(l0 == debugPoints(iDebugPoint)) THEN
1131 ! WRITE(*,*) 'Before Inviscid point ',l0
1132 ! WRITE(*,*) 'rho = ',ub(1), ' rhoT = ',gI1(1)
1133 ! WRITE(*,*) 'rhoU = ',ub(2),' rhoUT = ',gI1(2)
1134 ! WRITE(*,*) 'rhoV = ',ub(3),' rhoVT = ',gI1(3)
1135 ! WRITE(*,*) 'rhoE = ',ub(4),' rhoET = ',gI1(4)
1136 ! WRITE(*,*) 'U = ',bVelocity(1), ' uT = ',vTarget(1)
1137 ! WRITE(*,*) 'V = ',bVelocity(2), ' uT = ',vTarget(2)
1138 ! WRITE(*,*) 'uWall = ',vWall(1),' uWallTarget = ',vWallTarget(1)
1139 ! WRITE(*,*) 'vWall = ',vWall(2),' vWallTarget = ',vWallTarget(2)
1140 ! WRITE(*,*) 'PNORM = (',pnorm_vec,') vDotXi = ',vDotXi
1141 ! WRITE(*,*) 'normvec = ',norm_vec
1142 ! ENDIF
1143 ! ENDDO
1144 
1145  CALL sat_form_roe_matrices(nd, ub, gi1, gamma, pnorm_vec, tmat, tinv, lambda)
1146 
1147  specrad = 0.0_8
1148  ucon = vdotxi*xh
1149 
1150  ! ... only pick those Lambda that are incoming
1151  If ( sgn == 1 ) Then
1152  Do jj = 1, nd+2
1153  lambda(jj,jj) = max(lambda(jj,jj),0.0_8)
1154  specrad = max(abs(lambda(jj,jj)),abs(specrad))
1155  End Do
1156  Else
1157  Do jj = 1, nd+2
1158  lambda(jj,jj) = min(lambda(jj,jj),0.0_8)
1159  specrad = max(abs(lambda(jj,jj)),abs(specrad))
1160  End Do
1161  End If
1162 
1163 ! IF(TRIGGERED.eqv..FALSE.) THEN
1164 ! WRITE(*,*) 'XI_X_TILDE: ',XI_X_TILDE
1165 ! WRITE(*,*) 'XI_Y_TILDE: ',XI_Y_TILDE
1166 ! WRITE(*,*) 'XI_Z_TILDE: ',XI_Z_TILDE
1167 ! WRITE(*,*) 'XI_X: ',XI_X
1168 ! WRITE(*,*) 'XI_Y: ',XI_Y
1169 ! WRITE(*,*) 'XI_Z: ',XI_Z
1170 ! WRITE(*,*) 'BOUNDARYH: ',bndry_h
1171 ! WRITE(*,*) 'STATE: ',uB(:)
1172 ! WRITE(*,*) 'TARGET:',gI1(:)
1173 ! WRITE(*,*) 'UCON: ',uCon
1174 ! WRITE(*,*) 'SpecRad: ',specRad
1175 ! WRITE(*,*) 'vDotXi: ',vDotXi
1176 ! ENDIF
1177 ! debugPoint = 114273
1178 
1179  matx = matmul(lambda,tinv)
1180 
1181  ! ... compute the characteristic matrix
1182  aprime = matmul(tmat,matx)
1183 
1184  ! ... subtract off the target
1185  ub(1:nd+2) = ub(1:nd+2) - gi1(1:nd+2)
1186 
1187  ! ... compute the characteristic penalty vector
1188  gi2(1:nd+2) = matmul(aprime,ub)
1189 
1190  penaltyfac = sbpboundaryweight*sigmai1
1191  correction(:) = penaltyfac*gi2(:)
1192 
1193 ! DO iDebugPoint = 1, numDebugPoints
1194 ! IF(l0 == debugPoints(iDebugPoint)) THEN
1195 ! WRITE(*,*) 'Before Inviscid point ',l0
1196 ! WRITE(*,*) 'SIGMAI1 = ',sigmaI1,'bWeight = ',sbpBoundaryWeight
1197 ! WRITE(*,*) 'rhoRHS = ',rhoRHS(l0), ' dRho = ',correction(1)
1198 ! WRITE(*,*) 'rhoURHS = ',rhoVRHS(l0),' dRhoU = ',correction(2)
1199 ! WRITE(*,*) 'rhoVRHS = ',rhoVRHS(l0+numPointsBuffer),' dRhoV = ',correction(3)
1200 ! WRITE(*,*) 'rhoERHS = ',rhoERHS(l0),' dRhoE = ',correction(4)
1201 ! WRITE(*,*) 'cpRho = ',gI2(1)
1202 ! WRITE(*,*) 'cpRhoU = ',gI2(2)
1203 ! WRITE(*,*) 'cpRhoV = ',gI2(3)
1204 ! WRITE(*,*) 'cpRhoE = ',gI2(4)
1205 ! ENDIF
1206 ! ENDDO
1207 
1208  ! IBFAC would be considered here
1209  ! Inviscid correction
1210  rhorhs(l0) = rhorhs(l0) + ibfac_local * correction(1)
1211  DO idim = 1,numdim
1212  pointindex = l0 + (idim-1)*numpointsbuffer
1213  rhovrhs(pointindex) = rhovrhs(pointindex) + ibfac_local*correction(1+idim)
1214  END DO
1215  rhoerhs(l0) = rhoerhs(l0) + ibfac_local * correction(numdim+2)
1216 
1217 
1218 ! IF(TRIGGERED.eqv..FALSE.) THEN
1219 ! WRITE(*,*) 'CORRECTION(',l0,'): ',correction(:)
1220 ! TRIGGERED = .TRUE.
1221 ! ENDIF
1222  ! ... Species treatment (inviscid part)
1223  IF(numscalar > 0) THEN
1224  If (dsgn * ucon > 0) Then
1225  sfactor = -ucon * dsgn * ibfac_local * sbpboundaryweight
1226  Do jj = 1, numscalar
1227  pointindex = (jj-1)*numpointsbuffer + l0
1228  scalarrhs(pointindex) = scalarrhs(pointindex) + sfactor*scalarbuffer(pointindex)
1229  End do
1230  End If
1231  END IF
1232 
1233 
1234  ! VISCOUS PART
1235 
1236 ! uB(:) = cv(l0,:)
1237  ! MJA, repeated from above
1238  ub(1) = rhobuffer(l0)
1239  DO idim = 1, numdim
1240  ub(1+idim) = rhovbuffer(l0+(idim-1)*numpointsbuffer)
1241  END DO
1242  ub(numdim+2) = rhoebuffer(l0)
1243 
1244  ! ... wall temperature
1245  ! T_wall = bcic_WallTemp / ((gamref - 1.0_8)*tempref)
1246  ! new nondimensionalization...can be moved outside this loop altogether
1247  t_wall = bcwalltemp / tempref
1248 
1249  ! ... target state
1250  ! gI2(1) = RHO_TARGET
1251  gi2(1) = rhobuffer(l0)
1252 
1253  ! ... wall velocity
1254  ! gI2(2:ND+1) = XYZ_TAU(l0,1:ND) * RHO_TARGET
1255  ! gI2(2:ND+1) = 0.0_8
1256  ! MJA, non-moving meshes
1257  gi2(2:idim+1) = 0.0_8
1258 
1259  ! ... thermally perfect or not
1260  ! if (gas_dv_model == DV_MODEL_IDEALGAS) then
1261  ! gI2(ND+2) = RHO_TARGET * T_wall / GAMMA + 0.5_8 * SUM(gI2(2:ND+1)**2) / RHO_TARGET
1262 
1263  ! MJA, only implemented gas model right now
1264  ! old nondimen
1265  !gI2(numDim+2) = rhoTarget(l0) * T_wall / gamref + 0.5_8 * SUM(gI2(2:iDim+1)**2) / gI2(1)
1266  ! new dimen, roll spcGasConst and (gamma-1) into one const outside of loop
1267 
1268  ! use a temperature from the input file, if that temperature is <0,
1269  !use the target energy (temperature)
1270  IF (t_wall < 0) THEN
1271  gi2(numdim+2) = rhoetarget(l0)
1272  ELSE
1273  gi2(numdim+2) = gi2(1) * t_wall * spcgasconst / (gamref-1) &
1274  + 0.5_8 * sum(gi2(2:idim+1)**2) / gi2(1)
1275  ENDIF
1276 
1277  ! sigmaI2 = -(SAT_sigmaI2 / 4.0_8) * (SBP_invPdiag_block1(1)/bndry_h)**2 * &
1278  ! maxval(tv(l0,1:)) / RHO_TARGET * MAX(GAMMA / Pr, 5.0_8 / 3.0_8)
1279  ! sigmaI2 = sigmaI2 * REinv
1280  sigmai2 = -sat_sigmai2
1281 ! penaltyFac = (sbpBoundaryWeight/bndry_h)**2*sigmaI2/4.0_8 * &
1282 ! REInv*5.0_8/3.0_8*max(muBuffer(l0),lambdaBuffer(l0))/rhoTarget(l0)
1283 ! TESTING MTC
1284  penaltyfac = (sbpboundaryweight/bndry_h)**2*sigmai2/4.0_8 * &
1285  reinv*viscouspenaltyscale*max(mubuffer(l0),lambdabuffer(l0))/gi2(1)
1286 
1287 
1288  ! ... compute penalty
1289  ! ... compute penalty
1290  ! penalty(:) = sigmaI2 * (uB(1:ND+2) - gI2(1:ND+2))
1291  penalty(:) = penaltyfac * (ub(1:nd+2) - gi2(1:nd+2))
1292 
1293  ! ... add to the rhs
1294  ! do jj = 1, ND+2
1295  ! rhs(l0,jj) = rhs(l0,jj) + ibfac_local * penalty(jj)
1296  ! end do
1297  rhorhs(l0) = rhorhs(l0) + ibfac_local * penalty(1)
1298  DO idim = 1,numdim
1299  pointindex = l0 + (idim-1)*numpointsbuffer
1300  rhovrhs(pointindex) = rhovrhs(pointindex) + ibfac_local*penalty(1+idim)
1301  END DO
1302  rhoerhs(l0) = rhoerhs(l0) + ibfac_local * penalty(numdim+2)
1303 
1304 ! DO iDebugPoint = 1, numDebugPoints
1305 ! IF(l0 == debugPoints(iDebugPoint)) THEN
1306 ! WRITE(*,*) 'After Viscous point ',l0
1307 ! WRITE(*,*) 'rhoRHS = ',rhoRHS(l0), ' dRho = ',penalty(1)
1308 ! WRITE(*,*) 'rhoURHS = ',rhoVRHS(l0),' dRhoU = ',penalty(2)
1309 ! WRITE(*,*) 'rhoVRHS = ',rhoVRHS(l0+numPointsBuffer),' dRhoV = ',penalty(3)
1310 ! WRITE(*,*) 'rhoERHS = ',rhoERHS(l0),' dRhoE = ',penalty(4)
1311 ! WRITE(*,*) 'gI2(1) = ',gI2(1)
1312 ! WRITE(*,*) 'gI2(2) = ',gI2(2)
1313 ! WRITE(*,*) 'gI2(3) = ',gI2(3)
1314 ! WRITE(*,*) 'gI2(4) = ',gI2(4)
1315 ! ENDIF
1316 ! ENDDO
1317 
1318  !zero grad for species
1319  ! sigmaAux = ibfac_local * dsgn *1_8* SBP_invPdiag_block1(1) * JAC(l0)
1320  ! do jj = 1, nAuxVars
1321  ! ubAux = patch_ViscousFluxAux(jj,lp,normDir)
1322  ! rhs_auxVars(l0,jj) = rhs_auxVars(l0,jj) + sigmaAux* ubAux
1323  ! end do
1324 
1325  ! MJA
1326  ! species transport not fully implemented yet for viscous terms
1327 
1328  END DO
1329 
1330  DEALLOCATE(xix,vwall,bvelocity,vwalltarget,vtarget)
1331  DEALLOCATE(gi1, ub, tmat, tinv, penalty, correction)
1332  DEALLOCATE(gi2, aprime, lambda, matx, pnorm_vec)
1333 
1334  RETURN
1335 
1336  END SUBROUTINE noslip_isothermal
1337 
1338  SUBROUTINE roe_mat(numDim,inCV,outCV, &
1339  gamma,normDir,tMat,tInv,lambda,status)
1340 
1341  IMPLICIT NONE
1342 
1343  INTEGER(4), INTENT(IN) :: numdim
1344  REAL(8), INTENT(IN) :: incv(numdim+2), outcv(numdim+2)
1345  REAL(8), INTENT(IN) :: gamma, normdir(numdim)
1346  REAL(8), INTENT(OUT) :: tmat((numdim+2)*(numdim+2)),tinv((numdim+2)*(numdim+2)),lambda(numdim+2)
1347  INTEGER(4), INTENT(OUT) :: status
1348 
1349  REAL(8) :: gammam1,ucon,alpha,beta,phi2,theta,denomroe,xix,xiy,xiz
1350  REAL(8) :: xixtilde,xiytilde,xiztilde
1351  REAL(8) :: rhoin,rhov1in,rhov2in,rhov3in,rhom1in,rhoein
1352  REAL(8) :: rhoout,rhov1out,rhov2out,rhov3out,rhom1out,rhoeout
1353  REAL(8) :: inv1,inv2,inv3,outv1,outv2,outv3
1354  REAL(8) :: inke,inie,inp,inh
1355  REAL(8) :: outke,outie,outp,outh
1356  REAL(8) :: ccroe,bbroe,v1roe,v2roe,v3roe,hroe,keroe
1357  REAL(8) :: soundspeedroe
1358  INTEGER(4) :: idim
1359 
1360  status = 0
1361 
1362  IF(numdim == 1 .OR. numdim > 3) THEN
1363  status = 1
1364  RETURN
1365  ENDIF
1366 
1367  ! ... constants
1368  gammam1 = gamma - 1.0_8
1369 
1370  ! form primitives from "in"
1371  rhoin = incv(1)
1372  rhov1in = incv(2)
1373  rhov2in = incv(3)
1374  rhov3in = 0.0_8
1375  rhom1in = 1.0/rhoin
1376 
1377  inv1 = rhom1in*rhov1in
1378  inv2 = rhom1in*rhov2in
1379 
1380  IF(numdim == 2) THEN
1381  rhoein = incv(4)
1382  ELSE
1383  rhov3in = incv(4)
1384  rhoein = incv(5)
1385  ENDIF
1386  inv3 = rhom1in*rhov3in
1387 
1388  inke = 0.5_8 * (inv1**2 + inv2**2 + inv3**2)
1389  inie = rhom1in*rhoein - inke
1390  inp = gammam1*rhoin*inie
1391  inh = (gamma/gammam1)*inp/rhoin + inke
1392 
1393 
1394  ! form primitives from "out"
1395  rhoout = outcv(1)
1396  rhov1out = outcv(2)
1397  rhov2out = outcv(3)
1398  rhom1out = 1.0/rhoout
1399  outv1 = rhom1out*rhov1out
1400  outv2 = rhom1out*rhov2out
1401  IF(numdim == 2) THEN
1402  rhoeout = outcv(4)
1403  ELSE
1404  rhoeout = outcv(5)
1405  rhov3out = outcv(4)
1406  ENDIF
1407  outv3 = rhom1out*rhov3out
1408 
1409  outke = 0.5_8 * (outv1**2 + outv2**2 + outv3**2)
1410  outie = rhom1out*rhoeout - outke
1411  outp = gammam1*rhoout*outie
1412  outh = (gamma/gammam1)*outp/rhoout + outke
1413 
1414  ! Roe vars
1415  ccroe = sqrt(rhoout/rhoin)
1416  bbroe = 1.0_8 / (1.0_8 + ccroe)
1417  v1roe = (inv1 + outv1 * ccroe) * bbroe
1418  v2roe = (inv2 + outv2 * ccroe) * bbroe
1419  v3roe = (inv3 + outv3 * ccroe) * bbroe
1420  hroe = (inh + outh*ccroe)*bbroe
1421  keroe = 0.5_8 * (v1roe**2 + v2roe**2 + v3roe**2)
1422 
1423  IF((hroe - keroe) <= 0.0_8) THEN
1424  write (*,'(5(a,E20.8,1X))') "gamma = ", gamma, ", inKE = ", inke, ", inIE = ", inie, &
1425  ", inP = ", inp, ", inH = ", inh
1426  status = 1
1427  RETURN
1428  ENDIF
1429 
1430  soundspeedroe = sqrt(gammam1*(hroe - keroe))
1431 
1432  ! ... variable renaming
1433  xix = normdir(1)
1434  xiy = normdir(2)
1435  xiz = 0.0_8
1436  IF(numdim == 3) xiz = normdir(3)
1437  denomroe = 1.0_8 / sqrt(xix**2 + xiy**2 + xiz**2)
1438  xixtilde = xix * denomroe
1439  xiytilde = xiy * denomroe
1440  xiztilde = xiz * denomroe
1441 
1442  ucon = v1roe*xix + v2roe*xiy + v3roe*xiz
1443 
1444  ! ... directly from Pulliam & Chaussee
1445  alpha = 0.5_8*rhoin/soundspeedroe
1446  beta = 1.0_8/(rhoin*soundspeedroe)
1447  theta = xixtilde*v1roe + xiytilde*v2roe + xiztilde*v3roe
1448  phi2 = 0.5_8*(gammam1)*(v1roe**2 + v2roe**2 + v3roe**2)
1449 
1450  ! ... compute the diagonal matrix lambda
1451  lambda(:) = 0.0_8
1452  DO idim = 1, numdim
1453  lambda(idim) = ucon
1454  END DO
1455  lambda(numdim+1) = ucon + soundspeedroe * sqrt(xix**2 + xiy**2 + xiz**2)
1456  lambda(numdim+2) = ucon - soundspeedroe * sqrt(xix**2 + xiy**2 + xiz**2)
1457 
1458  if (numdim == 2) then
1459 
1460  tmat(1) = 1.0_8
1461  tmat(2) = v1roe
1462  tmat(3) = v2roe
1463  tmat(4) = phi2/gammam1
1464 
1465  tmat(5) = 0.0_8
1466  tmat(6) = xiytilde*rhoin
1467  tmat(7) = -xixtilde*rhoin
1468  tmat(8) = rhoin*(xiytilde*v1roe - xixtilde*v2roe)
1469 
1470  tmat(9) = alpha
1471  tmat(10) = alpha*(v1roe + xixtilde*soundspeedroe)
1472  tmat(11) = alpha*(v2roe + xiytilde*soundspeedroe)
1473  tmat(12) = alpha*((phi2 + soundspeedroe**2)/gammam1 + soundspeedroe*theta)
1474 
1475  tmat(13) = alpha
1476  tmat(14) = alpha*(v1roe - xixtilde*soundspeedroe)
1477  tmat(15) = alpha*(v2roe - xiytilde*soundspeedroe)
1478  tmat(16) = alpha*((phi2 + soundspeedroe**2)/gammam1 - soundspeedroe*theta)
1479 
1480  tinv(1) = 1.0_8 - phi2 / soundspeedroe**2
1481  tinv(2) = -(xiytilde * v1roe - xixtilde * v2roe) / rhoin
1482  tinv(3) = beta * (phi2 - soundspeedroe * theta)
1483  tinv(4) = beta * (phi2 + soundspeedroe * theta)
1484 
1485  tinv(5) = gammam1 * v1roe / soundspeedroe**2
1486  tinv(6) = xiytilde / rhoin
1487  tinv(7) = beta * (xixtilde * soundspeedroe - gammam1 * v1roe)
1488  tinv(8) = -beta * (xixtilde * soundspeedroe + gammam1 * v1roe)
1489 
1490  tinv(9) = gammam1 * v2roe / soundspeedroe**2
1491  tinv(10) = -xixtilde / rhoin
1492  tinv(11) = beta * (xiytilde * soundspeedroe - gammam1 * v2roe)
1493  tinv(12) = -beta * (xiytilde * soundspeedroe + gammam1 * v2roe)
1494 
1495  tinv(13) = -gammam1 / soundspeedroe**2
1496  tinv(14) = 0.0_8
1497  tinv(15) = beta * gammam1
1498  tinv(16) = beta * gammam1
1499 
1500  ELSE IF (numdim == 3) THEN
1501 
1502 
1503  tmat(1) = xixtilde
1504  tmat(2) = xixtilde * v1roe
1505  tmat(3) = xixtilde * v2roe + xiztilde * rhoin
1506  tmat(4) = xixtilde * v3roe - xiytilde * rhoin
1507  tmat(5) = xixtilde * phi2 / gammam1 + rhoin * (xiztilde * v2roe - xiytilde * v3roe)
1508 
1509  tmat(6) = xiytilde
1510  tmat(7) = xiytilde * v1roe - xiztilde * rhoin
1511  tmat(8) = xiytilde * v2roe
1512  tmat(9) = xiytilde * v3roe + xixtilde * rhoin
1513  tmat(10) = xiytilde * phi2 / gammam1 + rhoin * (xixtilde * v3roe - xiztilde * v1roe)
1514 
1515  tmat(11) = xiztilde
1516  tmat(12) = xiztilde * v1roe + xiytilde * rhoin
1517  tmat(13) = xiztilde * v2roe - xixtilde * rhoin
1518  tmat(14) = xiztilde * v3roe
1519  tmat(15) = xiztilde * phi2 / gammam1 + rhoin * (xiytilde * v1roe - xixtilde * v2roe)
1520 
1521  tmat(16) = alpha
1522  tmat(17) = alpha * (v1roe + xixtilde * soundspeedroe)
1523  tmat(18) = alpha * (v2roe + xiytilde * soundspeedroe)
1524  tmat(19) = alpha * (v3roe + xiztilde * soundspeedroe)
1525  tmat(20) = alpha * ((phi2 + soundspeedroe**2)/gammam1 + soundspeedroe * theta)
1526 
1527  tmat(21) = alpha
1528  tmat(22) = alpha * (v1roe - xixtilde * soundspeedroe)
1529  tmat(23) = alpha * (v2roe - xiytilde * soundspeedroe)
1530  tmat(24) = alpha * (v3roe - xiztilde * soundspeedroe)
1531  tmat(25) = alpha * ((phi2 + soundspeedroe**2)/gammam1 - soundspeedroe * theta)
1532 
1533  tinv(1) = xixtilde * (1.0_8 - phi2 / soundspeedroe**2) - (xiztilde * v2roe - xiytilde * v3roe) / rhoin
1534  tinv(2) = xiytilde * (1.0_8 - phi2 / soundspeedroe**2) - (xixtilde * v3roe - xiztilde * v1roe) / rhoin
1535  tinv(3) = xiztilde * (1.0_8 - phi2 / soundspeedroe**2) - (xiytilde * v1roe - xixtilde * v2roe) / rhoin
1536  tinv(4) = beta * (phi2 - soundspeedroe * theta)
1537  tinv(5) = beta * (phi2 + soundspeedroe * theta)
1538 
1539  tinv(6) = xixtilde * gammam1 * v1roe / soundspeedroe**2
1540  tinv(7) = -xiztilde / rhoin + xiytilde * gammam1 * v1roe / soundspeedroe**2
1541  tinv(8) = xiytilde / rhoin + xiztilde * gammam1 * v1roe / soundspeedroe**2
1542  tinv(9) = beta * (xixtilde * soundspeedroe - gammam1 * v1roe)
1543  tinv(10) = -beta * (xixtilde * soundspeedroe + gammam1 * v1roe)
1544 
1545  tinv(11) = xiztilde / rhoin + xixtilde * gammam1 * v2roe / soundspeedroe**2
1546  tinv(12) = xiytilde * gammam1 * v2roe / soundspeedroe**2
1547  tinv(13) = -xixtilde / rhoin + xiztilde * gammam1 * v2roe / soundspeedroe**2
1548  tinv(14) = beta * (xiytilde * soundspeedroe - gammam1 * v2roe)
1549  tinv(15) = -beta * (xiytilde * soundspeedroe + gammam1 * v2roe)
1550 
1551  tinv(16) = -xiytilde / rhoin + xixtilde * gammam1 * v3roe / soundspeedroe**2
1552  tinv(17) = xixtilde / rhoin + xiytilde * gammam1 * v3roe / soundspeedroe**2
1553  tinv(18) = xiztilde * gammam1 * v3roe / soundspeedroe**2
1554  tinv(19) = beta * (xiztilde * soundspeedroe - gammam1 * v3roe)
1555  tinv(20) = -beta * (xiztilde * soundspeedroe + gammam1 * v3roe)
1556 
1557  tinv(21) = -xixtilde * gammam1 / soundspeedroe**2
1558  tinv(22) = -xiytilde * gammam1 / soundspeedroe**2
1559  tinv(23) = -xiztilde * gammam1 / soundspeedroe**2
1560  tinv(24) = beta * gammam1
1561  tinv(25) = beta * gammam1
1562 
1563 
1564  END IF
1565 
1566  END SUBROUTINE roe_mat
1567 
1568 
1569 
1570 
1571  subroutine sat_form_roe_matrices_tp(u_in, u_out, ND, gamma, gamref, norm, tmat, tinv, lambda)
1572 
1573 ! USE ModGlobal
1574 ! USE ModMPI
1575  IMPLICIT NONE
1576 
1577  Real(8), Pointer :: u_in(:), u_out(:), tmat(:,:), tinv(:,:), norm(:), lambda(:,:)
1578  Real(8) :: gamma, rho_in, ux_in, uy_in, uz_in, ke_in, p_in, en_in, h_in
1579  Real(8) :: gm1, rho_out, ux_out, uy_out, uz_out, ke_out, p_out, en_out, h_out
1580  Real(8) :: cc, bb, ux_roe, uy_roe, uz_roe, h_roe, ke_roe, spd_snd_roe, t_roe, en_roe
1581  Real(8) :: xi_x_tilde, xi_y_tilde, xi_z_tilde, ucon, rho, alpha, theta, beta, phi2
1582  Real(8) :: xi_x, xi_y, xi_z, xi_t, denom, gamref, up, um
1583 
1584  Integer :: n, nd, jj, i, j, counter,ii
1585  Logical, dimension(ND+2) :: bwork
1586  Integer, dimension(ND+2) :: order
1587 
1588  ! ... Lapack arguments
1589  CHARACTER(len=1) :: jobvl, jobvr
1590  INTEGER :: info, lda, ldvl, ldvr, lwork,m
1591  Real(8), dimension(ND+2,ND+2) :: a, vl, vr,b,lam, pmat
1592  Real(8), dimension(ND+2) :: wi, wr
1593  Real(8), dimension(300) :: work
1594  ! ... for right vector inversion
1595  integer, dimension(ND+2) :: ipiv
1596 
1597  gamma = gamref
1598 
1599  ! ... temporary
1600  xi_t = 0.0_8
1601 
1602  jobvl = 'N'
1603  jobvr = 'V'
1604  lda = nd + 2
1605  ldvl = nd + 2
1606  ldvr = nd + 2
1607  m = nd+2
1608 
1609  lwork = 300
1610 
1611  lam(:,:) = 0.0_8
1612  a(:,:) = 0.0_8
1613  vl(:,:) = 0.0_8
1614  vr(:,:) = 0.0_8
1615  wr(:) = 0.0_8
1616  wi(:) = 0.0_8
1617  info = 0
1618 
1619 
1620  ! ... constants
1621  gm1 = gamma - 1.0_8
1622 
1623  ! ... form primative variables from "in"
1624  rho_in = u_in(1)
1625  ux_in = u_in(2) / rho_in; uy_in = 0.0_8; uz_in = 0.0_8
1626  if (nd >= 2) uy_in = u_in(3) / rho_in
1627  if (nd == 3) uz_in = u_in(4) / rho_in
1628  ke_in = 0.5_8 * (ux_in**2 + uy_in**2 + uz_in**2)
1629  en_in = u_in(nd+2)/rho_in - ke_in
1630  p_in = gm1 * rho_in * en_in
1631  h_in = (gamma / gm1) * p_in / rho_in + ke_in
1632 
1633  ! ... form primative variables from "out"
1634  rho_out = u_out(1)
1635  ux_out = u_out(2) / rho_out; uy_out = 0.0_8; uz_out = 0.0_8
1636  if (nd >= 2) uy_out = u_out(3) / rho_out
1637  if (nd == 3) uz_out = u_out(4) / rho_out
1638  ke_out = 0.5_8 * (ux_out**2 + uy_out**2 + uz_out**2)
1639  en_out = u_out(nd+2)/rho_out - ke_out
1640  p_out = gm1 * rho_out * en_out
1641  h_out = (gamma / gm1) * p_out / rho_out + ke_out
1642 
1643  ! ... form Roe variables
1644  cc = sqrt(rho_out / rho_in)
1645  bb = 1.0_8 / (1.0_8 + cc)
1646  ux_roe = (ux_in + ux_out * cc) * bb
1647  uy_roe = (uy_in + uy_out * cc) * bb
1648  uz_roe = (uz_in + uz_out * cc) * bb
1649  h_roe = ( h_in + h_out * cc) * bb
1650  ke_roe = 0.5_8 * (ux_roe**2 + uy_roe**2 + uz_roe**2)
1651  spd_snd_roe = sqrt(gamma/gamref * gm1*(h_roe - ke_roe))
1652 
1653  ! ... temporarily
1654  t_roe = gamma*(en_in + en_out * cc) * bb
1655  en_roe = t_roe/gamma
1656 
1657  ! ... variable renaming
1658  xi_x = norm(1); xi_y = 0.0_8; xi_z = 0.0_8
1659  if (nd >= 2) xi_y = norm(2)
1660  if (nd == 3) xi_z = norm(3)
1661  denom = 1.0_8 / sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1662  xi_x_tilde = xi_x * denom
1663  xi_y_tilde = xi_y * denom
1664  xi_z_tilde = xi_z * denom
1665  rho = rho_in
1666 
1667  ! ... directly from Pulliam & Chaussee
1668  theta = xi_x * ux_roe + xi_y * uy_roe + xi_z * uz_roe
1669  phi2 = 0.5_8 * (gamma - 1.0_8) * (ux_roe**2 + uy_roe**2 + uz_roe**2)
1670 
1671  if (nd == 2) then
1672 
1673  a(1,1) = xi_t
1674  a(1,2) = xi_x
1675  a(1,3) = xi_y
1676  a(1,4) = 0.0_8
1677 
1678  a(2,1) = xi_x * ( t_roe / gamref * (gamref - 1.0_8) - (gm1 * en_roe - phi2)) - theta * ux_roe
1679  a(2,2) = xi_t + theta - xi_x * (gamma - 2.0_8) * ux_roe
1680  a(2,3) = xi_y * ux_roe - gm1 * xi_x * uy_roe
1681  a(2,4) = gm1 * xi_x
1682 
1683  a(3,1) = xi_y * ( t_roe / gamref * (gamref - 1.0_8) - (gm1 * en_roe - phi2)) - theta * uy_roe
1684  a(3,2) = xi_x * uy_roe - gm1 * xi_y * ux_roe
1685  a(3,3) = xi_t + theta - xi_y * (gamma - 2.0_8) * uy_roe
1686  a(3,4) = gm1 * xi_y
1687 
1688  a(4,1) = theta * (2.0_8 * phi2 - gamma * (en_roe + ke_roe))
1689  a(4,2) = xi_x * (t_roe / gamref * (gamref - 1.0_8) + en_roe + phi2 / gm1 ) - gm1 * theta * ux_roe
1690  a(4,3) = xi_y * (t_roe / gamref * (gamref - 1.0_8) + en_roe + phi2 / gm1 ) - gm1 * theta * uy_roe
1691  a(4,4) = xi_t + gamma * theta
1692 
1693  else if (nd == 3) then
1694 
1695  a(1,1) = xi_t
1696  a(1,2) = xi_x
1697  a(1,3) = xi_y
1698  a(1,4) = xi_z
1699  a(1,5) = 0.0_8
1700 
1701  a(2,1) = xi_x * ( t_roe / gamref * (gamref - 1.0_8) - (gm1 * en_roe - phi2)) - theta * ux_roe
1702  a(2,2) = xi_t + theta - xi_x * (gamma - 2.0_8) * ux_roe
1703  a(2,3) = xi_y * ux_roe - gm1 * xi_x * uy_roe
1704  a(2,4) = xi_z * ux_roe - gm1 * xi_x * uz_roe
1705  a(2,5) = gm1 * xi_x
1706 
1707  a(3,1) = xi_y * ( t_roe / gamref * (gamref - 1.0_8) - (gm1 * en_roe - phi2)) - theta * uy_roe
1708  a(3,2) = xi_x * uy_roe - gm1 * xi_y * ux_roe
1709  a(3,3) = xi_t + theta - xi_y * (gamma - 2.0_8) * uy_roe
1710  a(3,4) = xi_z * uy_roe - gm1 * xi_y * uz_roe
1711  a(3,5) = gm1 * xi_y
1712 
1713  a(4,1) = xi_z * ( t_roe / gamref * (gamref - 1.0_8) - (gm1 * en_roe - phi2)) - theta * uz_roe
1714  a(4,2) = xi_x * uz_roe - gm1 * xi_z * ux_roe
1715  a(4,3) = xi_y * uz_roe - gm1 * xi_z * uy_roe
1716  a(4,4) = xi_t + theta - xi_z * (gamma - 2.0_8) * uz_roe
1717  a(4,5) = gm1 * xi_z
1718 
1719  a(5,1) = theta * (2.0_8 * phi2 - gamma * (en_roe + ke_roe))
1720  a(5,2) = xi_x * (t_roe / gamref * (gamref - 1.0_8) + en_roe + phi2 / gm1 ) - gm1 * theta * ux_roe
1721  a(5,3) = xi_y * (t_roe / gamref * (gamref - 1.0_8) + en_roe + phi2 / gm1 ) - gm1 * theta * uy_roe
1722  a(5,4) = xi_z * (t_roe / gamref * (gamref - 1.0_8) + en_roe + phi2 / gm1 ) - gm1 * theta * uz_roe
1723  a(5,5) = xi_t + gamma * theta
1724 
1725  end if
1726 
1727  ! ... find the Eigenvalues and Eigenvectors
1728  ! ... Lapack driver DGEEV
1729 !#ifdef HAVE_BLAS
1730 ! CALL DGEEV( JOBVL, JOBVR, M, A, LDA, WR, WI, VL, LDVL, VR, &
1731 ! LDVR, WORK, LWORK, INFO )
1732 !#endif
1733 
1734  lambda(:,:) = 0.0_8
1735  pmat(:,:) = 0.0_8
1736 
1737  ! ... form transformation matrices
1738  ucon = ux_roe * xi_x + uy_roe * xi_y + uz_roe * xi_z
1739  up = ucon + 0.1_8 * spd_snd_roe * sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1740  um = ucon - 0.1_8 * spd_snd_roe * sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1741  counter = 0
1742 
1743  ! ... Sort eigenvalues: {U,U,U,U+c,U-c}
1744  do jj = 1, nd+2
1745 
1746  if(wr(jj) .lt. up .and. wr(jj) .gt. um) then
1747 
1748  counter = counter + 1
1749  lambda(counter,counter) = wr(jj)
1750  order(jj) = counter
1751  pmat(jj,counter) = 1.0_8
1752  elseif(wr(jj) .gt. ucon) then
1753  lambda(nd+1,nd+1) = wr(jj)
1754  order(jj) = nd+1
1755  pmat(jj,nd+1) = 1.0_8
1756  elseif(wr(jj) .lt. ucon) then
1757  lambda(nd+2,nd+2) = wr(jj)
1758  order(jj) = nd+2
1759  pmat(jj,nd+2) = 1.0_8
1760  end if
1761  end do
1762 
1763  tmat = matmul(vr,pmat)
1764  vr = tmat
1765 
1766  ! ... invert VR
1767  ! ... using a pseudoinvert using singular value decomposition
1768  ! ... VR = U*D*V^T
1769 
1770  ! ... 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
1771  jobvl = 'A'
1772  wi(:) = 0.0_8
1773  vl(:,:) = 0.0_8
1774  a(:,:) = 0.0_8
1775 
1776  ! ... VR = VL*diag(WI)*A^T
1777 !#ifdef HAVE_BLAS
1778 ! call DGESVD(JOBVL,JOBVL,LDVR,LDVR,VR,LDVR,WI,VL,LDVR,A,LDVR,WORK,LWORK,INFO)
1779 !#endif
1780 
1781  tinv(:,:) = 0.0_8
1782  ! ... 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
1783  do i = 1,nd+2
1784  if(wi(i) .gt. tiny(wi(i))) then
1785  tinv(i,i) = 1.0_8/wi(i)
1786  else
1787  tinv(i,i) = 0.0_8
1788  end if
1789  end do
1790 
1791 
1792  tinv = matmul(transpose(a),tinv)
1793  tinv = matmul(tinv,transpose(vl))
1794 
1795  end subroutine sat_form_roe_matrices_tp
1796 
1797 
1798 
1799  subroutine sat_form_roe_matrices(ND, u_in, u_out, gamma, norm, tmat, tinv, lambda)
1800 
1801 ! USE ModGlobal
1802 
1803  IMPLICIT NONE
1804 
1805  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)
1806  Real(8) :: gamma, rho_in, ux_in, uy_in, uz_in, ke_in, p_in, en_in, h_in
1807  Real(8) :: gm1, rho_out, ux_out, uy_out, uz_out, ke_out, p_out, en_out, h_out
1808  Real(8) :: cc, bb, ux_roe, uy_roe, uz_roe, h_roe, ke_roe, spd_snd_roe
1809  Real(8) :: xi_x_tilde, xi_y_tilde, xi_z_tilde, ucon, rho, alpha, theta, beta, phi2
1810  Real(8) :: xi_x, xi_y, xi_z, denom
1811  Integer :: n, nd, jj
1812 
1813  ! ... size of this problem
1814  !ND = size(u_in) - 2
1815 
1816  ! ... constants
1817  gm1 = gamma - 1.0_8
1818 
1819  ! ... form primative variables from "in"
1820  rho_in = u_in(1)
1821  ux_in = u_in(2) / rho_in; uy_in = 0.0_8; uz_in = 0.0_8
1822  if (nd >= 2) uy_in = u_in(3) / rho_in
1823  if (nd == 3) uz_in = u_in(4) / rho_in
1824  ke_in = 0.5_8 * (ux_in**2 + uy_in**2 + uz_in**2)
1825  en_in = u_in(nd+2)/rho_in - ke_in
1826  p_in = gm1 * rho_in * en_in
1827  h_in = (gamma / gm1) * p_in / rho_in + ke_in
1828 
1829  ! ... form primative variables from "out"
1830  rho_out = u_out(1)
1831  ux_out = u_out(2) / rho_out; uy_out = 0.0_8; uz_out = 0.0_8
1832  if (nd >= 2) uy_out = u_out(3) / rho_out
1833  if (nd == 3) uz_out = u_out(4) / rho_out
1834  ke_out = 0.5_8 * (ux_out**2 + uy_out**2 + uz_out**2)
1835  en_out = u_out(nd+2)/rho_out - ke_out
1836  p_out = gm1 * rho_out * en_out
1837  h_out = (gamma / gm1) * p_out / rho_out + ke_out
1838 
1839  ! ... form Roe variables
1840  cc = sqrt(rho_out / rho_in)
1841  bb = 1.0_8 / (1.0_8 + cc)
1842  ux_roe = (ux_in + ux_out * cc) * bb
1843  uy_roe = (uy_in + uy_out * cc) * bb
1844  uz_roe = (uz_in + uz_out * cc) * bb
1845  h_roe = ( h_in + h_out * cc) * bb
1846  ke_roe = 0.5_8 * (ux_roe**2 + uy_roe**2 + uz_roe**2)
1847  if (h_roe - ke_roe <= 0.0_8) then
1848  ! write (*,'(6(E20.8,1X))') gamma, h_Roe, ke_Roe, ux_roe, uy_roe, uz_roe
1849  write (*,'(5(a,E20.8,1X))') "SATUTIL:FORM_ROE_MATRICES: ERRONEOUS CONDITIONS: gamma = ", gamma, ", ke_in = ", &
1850  ke_in, ", en_in = ", en_in, ", p_in = ", p_in, ", h_in = ", h_in
1851  stop
1852  end if
1853  spd_snd_roe = sqrt(gm1*(h_roe - ke_roe))
1854 
1855  ! ... variable renaming
1856  xi_x = norm(1); xi_y = 0.0_8; xi_z = 0.0_8
1857  if (nd >= 2) xi_y = norm(2)
1858  if (nd == 3) xi_z = norm(3)
1859  denom = 1.0_8 / sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1860  xi_x_tilde = xi_x * denom
1861  xi_y_tilde = xi_y * denom
1862  xi_z_tilde = xi_z * denom
1863  rho = rho_in
1864 
1865  ! ... form transformation matrices
1866  ucon = ux_roe * xi_x + uy_roe * xi_y + uz_roe * xi_z
1867 
1868  if (nd == 2) then
1869 
1870  ! ... directly from Pulliam & Chaussee
1871  alpha = 0.5_8 * rho / spd_snd_roe
1872  beta = 1.0_8 / (rho * spd_snd_roe)
1873  theta = xi_x_tilde * ux_roe + xi_y_tilde * uy_roe
1874  phi2 = 0.5_8 * (gamma - 1.0_8) * (ux_roe**2 + uy_roe**2)
1875  tmat(1,1) = 1.0_8
1876  tmat(1,2) = 0.0_8
1877  tmat(1,3) = alpha
1878  tmat(1,4) = alpha
1879  tmat(2,1) = ux_roe
1880  tmat(2,2) = xi_y_tilde * rho
1881  tmat(2,3) = alpha * (ux_roe + xi_x_tilde * spd_snd_roe)
1882  tmat(2,4) = alpha * (ux_roe - xi_x_tilde * spd_snd_roe)
1883  tmat(3,1) = uy_roe
1884  tmat(3,2) = -xi_x_tilde * rho
1885  tmat(3,3) = alpha * (uy_roe + xi_y_tilde * spd_snd_roe)
1886  tmat(3,4) = alpha * (uy_roe - xi_y_tilde * spd_snd_roe)
1887  tmat(4,1) = phi2 / (gamma - 1.0_8)
1888  tmat(4,2) = rho * (xi_y_tilde * ux_roe - xi_x_tilde * uy_roe)
1889  tmat(4,3) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) + spd_snd_roe * theta)
1890  tmat(4,4) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) - spd_snd_roe * theta)
1891 
1892  tinv(1,1) = 1.0_8 - phi2 / spd_snd_roe**2
1893  tinv(1,2) = (gamma - 1.0_8) * ux_roe / spd_snd_roe**2
1894  tinv(1,3) = (gamma - 1.0_8) * uy_roe / spd_snd_roe**2
1895  tinv(1,4) = -(gamma - 1.0_8) / spd_snd_roe**2
1896  tinv(2,1) = -(xi_y_tilde * ux_roe - xi_x_tilde * uy_roe) / rho
1897  tinv(2,2) = xi_y_tilde / rho
1898  tinv(2,3) = -xi_x_tilde / rho
1899  tinv(2,4) = 0.0_8
1900  tinv(3,1) = beta * (phi2 - spd_snd_roe * theta)
1901  tinv(3,2) = beta * (xi_x_tilde * spd_snd_roe - (gamma - 1.0_8) * ux_roe)
1902  tinv(3,3) = beta * (xi_y_tilde * spd_snd_roe - (gamma - 1.0_8) * uy_roe)
1903  tinv(3,4) = beta * (gamma - 1.0_8)
1904  tinv(4,1) = beta * (phi2 + spd_snd_roe * theta)
1905  tinv(4,2) = -beta * (xi_x_tilde * spd_snd_roe + (gamma - 1.0_8) * ux_roe)
1906  tinv(4,3) = -beta * (xi_y_tilde * spd_snd_roe + (gamma - 1.0_8) * uy_roe)
1907  tinv(4,4) = beta * (gamma - 1.0_8)
1908 
1909  ! ... compute the diagonal matrix lambda
1910  lambda(:,:) = 0.0_8
1911  do jj = 1, nd
1912  lambda(jj,jj) = ucon
1913  end do
1914  lambda(nd+1,nd+1) = ucon + spd_snd_roe * sqrt(xi_x**2 + xi_y**2)
1915  lambda(nd+2,nd+2) = ucon - spd_snd_roe * sqrt(xi_x**2 + xi_y**2)
1916 
1917  else if (nd == 3) then
1918 
1919  ! ... directly from Pulliam & Chaussee
1920  alpha = 0.5_8 * rho / spd_snd_roe
1921  beta = 1.0_8 / (rho * spd_snd_roe)
1922  theta = xi_x_tilde * ux_roe + xi_y_tilde * uy_roe + xi_z_tilde * uz_roe
1923  phi2 = 0.5_8 * (gamma - 1.0_8) * (ux_roe**2 + uy_roe**2 + uz_roe**2)
1924 
1925  tmat(1,1) = xi_x_tilde
1926  tmat(1,2) = xi_y_tilde
1927  tmat(1,3) = xi_z_tilde
1928  tmat(1,4) = alpha
1929  tmat(1,5) = alpha
1930  tmat(2,1) = xi_x_tilde * ux_roe
1931  tmat(2,2) = xi_y_tilde * ux_roe - xi_z_tilde * rho
1932  tmat(2,3) = xi_z_tilde * ux_roe + xi_y_tilde * rho
1933  tmat(2,4) = alpha * (ux_roe + xi_x_tilde * spd_snd_roe)
1934  tmat(2,5) = alpha * (ux_roe - xi_x_tilde * spd_snd_roe)
1935  tmat(3,1) = xi_x_tilde * uy_roe + xi_z_tilde * rho
1936  tmat(3,2) = xi_y_tilde * uy_roe
1937  tmat(3,3) = xi_z_tilde * uy_roe - xi_x_tilde * rho
1938  tmat(3,4) = alpha * (uy_roe + xi_y_tilde * spd_snd_roe)
1939  tmat(3,5) = alpha * (uy_roe - xi_y_tilde * spd_snd_roe)
1940  tmat(4,1) = xi_x_tilde * uz_roe - xi_y_tilde * rho
1941  tmat(4,2) = xi_y_tilde * uz_roe + xi_x_tilde * rho
1942  tmat(4,3) = xi_z_tilde * uz_roe
1943  tmat(4,4) = alpha * (uz_roe + xi_z_tilde * spd_snd_roe)
1944  tmat(4,5) = alpha * (uz_roe - xi_z_tilde * spd_snd_roe)
1945  tmat(5,1) = xi_x_tilde * phi2 / (gamma - 1.0_8) + rho * (xi_z_tilde * uy_roe - xi_y_tilde * uz_roe)
1946  tmat(5,2) = xi_y_tilde * phi2 / (gamma - 1.0_8) + rho * (xi_x_tilde * uz_roe - xi_z_tilde * ux_roe)
1947  tmat(5,3) = xi_z_tilde * phi2 / (gamma - 1.0_8) + rho * (xi_y_tilde * ux_roe - xi_x_tilde * uy_roe)
1948  tmat(5,4) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) + spd_snd_roe * theta)
1949  tmat(5,5) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) - spd_snd_roe * theta)
1950 
1951  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
1952  tinv(1,2) = xi_x_tilde * (gamma - 1.0_8) * ux_roe / spd_snd_roe**2
1953  tinv(1,3) = xi_z_tilde / rho + xi_x_tilde * (gamma - 1.0_8) * uy_roe / spd_snd_roe**2
1954  tinv(1,4) = -xi_y_tilde / rho + xi_x_tilde * (gamma - 1.0_8) * uz_roe / spd_snd_roe**2
1955  tinv(1,5) = -xi_x_tilde * (gamma - 1.0_8) / spd_snd_roe**2
1956  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
1957  tinv(2,2) = -xi_z_tilde / rho + xi_y_tilde * (gamma - 1.0_8) * ux_roe / spd_snd_roe**2
1958  tinv(2,3) = xi_y_tilde * (gamma - 1.0_8) * uy_roe / spd_snd_roe**2
1959  tinv(2,4) = xi_x_tilde / rho + xi_y_tilde * (gamma - 1.0_8) * uz_roe / spd_snd_roe**2
1960  tinv(2,5) = -xi_y_tilde * (gamma - 1.0_8) / spd_snd_roe**2
1961  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
1962  tinv(3,2) = xi_y_tilde / rho + xi_z_tilde * (gamma - 1.0_8) * ux_roe / spd_snd_roe**2
1963  tinv(3,3) = -xi_x_tilde / rho + xi_z_tilde * (gamma - 1.0_8) * uy_roe / spd_snd_roe**2
1964  tinv(3,4) = xi_z_tilde * (gamma - 1.0_8) * uz_roe / spd_snd_roe**2
1965  tinv(3,5) = -xi_z_tilde * (gamma - 1.0_8) / spd_snd_roe**2
1966  tinv(4,1) = beta * (phi2 - spd_snd_roe * theta)
1967  tinv(4,2) = beta * (xi_x_tilde * spd_snd_roe - (gamma - 1.0_8) * ux_roe)
1968  tinv(4,3) = beta * (xi_y_tilde * spd_snd_roe - (gamma - 1.0_8) * uy_roe)
1969  tinv(4,4) = beta * (xi_z_tilde * spd_snd_roe - (gamma - 1.0_8) * uz_roe)
1970  tinv(4,5) = beta * (gamma - 1.0_8)
1971  tinv(5,1) = beta * (phi2 + spd_snd_roe * theta)
1972  tinv(5,2) = -beta * (xi_x_tilde * spd_snd_roe + (gamma - 1.0_8) * ux_roe)
1973  tinv(5,3) = -beta * (xi_y_tilde * spd_snd_roe + (gamma - 1.0_8) * uy_roe)
1974  tinv(5,4) = -beta * (xi_z_tilde * spd_snd_roe + (gamma - 1.0_8) * uz_roe)
1975  tinv(5,5) = beta * (gamma - 1.0_8)
1976 
1977  ! ... compute the diagonal matrix lambda
1978  lambda(:,:) = 0.0_8
1979  do jj = 1, nd
1980  lambda(jj,jj) = ucon
1981  end do
1982  lambda(nd+1,nd+1) = ucon + spd_snd_roe * sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1983  lambda(nd+2,nd+2) = ucon - spd_snd_roe * sqrt(xi_x**2 + xi_y**2 + xi_z**2)
1984 
1985  end if
1986 
1987  end subroutine sat_form_roe_matrices
1988 
1989  subroutine sat_form_roe_matrices_test(ND, u_in, u_out, gamma, norm, metricData,tmat, tinv, lambda)
1990 
1991  ! USE ModGlobal
1992 
1993  IMPLICIT NONE
1994 
1995  Real(8) :: u_in(nd+2), u_out(nd+2), tmat(nd+2,nd+2), norm(nd)
1996  REAL(8) :: metricData(9), tinv(nd+2,nd+2), lambda(nd+2,nd+2)
1997  Real(8) :: gamma, rho_in, ux_in, uy_in, uz_in, ke_in, p_in, en_in, h_in
1998  Real(8) :: gm1, rho_out, ux_out, uy_out, uz_out, ke_out, p_out, en_out, h_out
1999  Real(8) :: cc, bb, ux_Roe, uy_Roe, uz_Roe, h_Roe, ke_Roe, spd_snd_Roe
2000  Real(8) :: XI_X_TILDE, XI_Y_TILDE, XI_Z_TILDE, ucon, rho, alpha, theta, beta, phi2
2001  Real(8) :: XI_X, XI_Y, XI_Z, b1, b2, U2, denom
2002  Integer :: N, ND, jj, dir, i
2003  Real(8) :: uhat, vhat, what, k_x, k_y, k_z, l_x, l_y, l_z, m_x, m_y, m_z, ja!,check(5,5)
2004 
2005  ! ... size of this problem
2006  !ND = size(u_in) - 2
2007 
2008  ! ... constants
2009  gm1 = gamma - 1.0_8
2010 
2011  ! ... form primative variables from "in"
2012  rho_in = u_in(1)
2013  ux_in = u_in(2) / rho_in; uy_in = 0.0_8; uz_in = 0.0_8
2014  if (nd >= 2) uy_in = u_in(3) / rho_in
2015  if (nd == 3) uz_in = u_in(4) / rho_in
2016  ke_in = 0.5_8 * (ux_in**2 + uy_in**2 + uz_in**2)
2017  en_in = u_in(nd+2)/rho_in - ke_in
2018  p_in = gm1 * rho_in * en_in
2019  h_in = (gamma / gm1) * p_in / rho_in + ke_in
2020 
2021  ! ... form primative variables from "out"
2022  rho_out = u_out(1)
2023  ux_out = u_out(2) / rho_out; uy_out = 0.0_8; uz_out = 0.0_8
2024  if (nd >= 2) uy_out = u_out(3) / rho_out
2025  if (nd == 3) uz_out = u_out(4) / rho_out
2026  ke_out = 0.5_8 * (ux_out**2 + uy_out**2 + uz_out**2)
2027  en_out = u_out(nd+2)/rho_out - ke_out
2028  p_out = gm1 * rho_out * en_out
2029  h_out = (gamma / gm1) * p_out / rho_out + ke_out
2030 
2031  ! ... form Roe variables
2032  cc = sqrt(rho_out / rho_in)
2033  bb = 1.0_8 / (1.0_8 + cc)
2034  ux_roe = (ux_in + ux_out * cc) * bb
2035  uy_roe = (uy_in + uy_out * cc) * bb
2036  uz_roe = (uz_in + uz_out * cc) * bb
2037  h_roe = ( h_in + h_out * cc) * bb
2038  ke_roe = 0.5_8 * (ux_roe**2 + uy_roe**2 + uz_roe**2)
2039  if (h_roe - ke_roe <= 0.0_8) then
2040  ! write (*,'(6(E20.8,1X))') gamma, h_Roe, ke_Roe, ux_roe, uy_roe, uz_roe
2041  write (*,'(5(a,E20.8,1X))') "SATUTIL:FORM_ROE_MATRICES: ERRONEOUS CONDITIONS: gamma = ", gamma, ", ke_in = ", &
2042  ke_in, ", en_in = ", en_in, ", p_in = ", p_in, ", h_in = ", h_in
2043  stop
2044  end if
2045  spd_snd_roe = sqrt(gm1*(h_roe - ke_roe))
2046 
2047  ! ... variable renaming
2048  xi_x = norm(1); xi_y = 0.0_8; xi_z = 0.0_8
2049  if (nd >= 2) xi_y = norm(2)
2050  if (nd == 3) xi_z = norm(3)
2051  denom = 1.0_8 / sqrt(xi_x**2 + xi_y**2 + xi_z**2)
2052  xi_x_tilde = xi_x * denom
2053  xi_y_tilde = xi_y * denom
2054  xi_z_tilde = xi_z * denom
2055  rho = rho_in
2056 
2057  ! ... form transformation matrices
2058  ucon = ux_roe * xi_x + uy_roe * xi_y + uz_roe * xi_z
2059 
2060  if (nd == 2) then
2061 
2062  ! ... directly from Pulliam & Chaussee
2063  alpha = 0.5_8 * rho / spd_snd_roe
2064  beta = 1.0_8 / (rho * spd_snd_roe)
2065  theta = xi_x_tilde * ux_roe + xi_y_tilde * uy_roe
2066  phi2 = 0.5_8 * (gamma - 1.0_8) * (ux_roe**2 + uy_roe**2)
2067  tmat(1,1) = 1.0_8
2068  tmat(1,2) = 0.0_8
2069  tmat(1,3) = alpha
2070  tmat(1,4) = alpha
2071  tmat(2,1) = ux_roe
2072  tmat(2,2) = xi_y_tilde * rho
2073  tmat(2,3) = alpha * (ux_roe + xi_x_tilde * spd_snd_roe)
2074  tmat(2,4) = alpha * (ux_roe - xi_x_tilde * spd_snd_roe)
2075  tmat(3,1) = uy_roe
2076  tmat(3,2) = -xi_x_tilde * rho
2077  tmat(3,3) = alpha * (uy_roe + xi_y_tilde * spd_snd_roe)
2078  tmat(3,4) = alpha * (uy_roe - xi_y_tilde * spd_snd_roe)
2079  tmat(4,1) = phi2 / (gamma - 1.0_8)
2080  tmat(4,2) = rho * (xi_y_tilde * ux_roe - xi_x_tilde * uy_roe)
2081  tmat(4,3) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) + spd_snd_roe * theta)
2082  tmat(4,4) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) - spd_snd_roe * theta)
2083 
2084  tinv(1,1) = 1.0_8 - phi2 / spd_snd_roe**2
2085  tinv(1,2) = (gamma - 1.0_8) * ux_roe / spd_snd_roe**2
2086  tinv(1,3) = (gamma - 1.0_8) * uy_roe / spd_snd_roe**2
2087  tinv(1,4) = -(gamma - 1.0_8) / spd_snd_roe**2
2088  tinv(2,1) = -(xi_y_tilde * ux_roe - xi_x_tilde * uy_roe) / rho
2089  tinv(2,2) = xi_y_tilde / rho
2090  tinv(2,3) = -xi_x_tilde / rho
2091  tinv(2,4) = 0.0_8
2092  tinv(3,1) = beta * (phi2 - spd_snd_roe * theta)
2093  tinv(3,2) = beta * (xi_x_tilde * spd_snd_roe - (gamma - 1.0_8) * ux_roe)
2094  tinv(3,3) = beta * (xi_y_tilde * spd_snd_roe - (gamma - 1.0_8) * uy_roe)
2095  tinv(3,4) = beta * (gamma - 1.0_8)
2096  tinv(4,1) = beta * (phi2 + spd_snd_roe * theta)
2097  tinv(4,2) = -beta * (xi_x_tilde * spd_snd_roe + (gamma - 1.0_8) * ux_roe)
2098  tinv(4,3) = -beta * (xi_y_tilde * spd_snd_roe + (gamma - 1.0_8) * uy_roe)
2099  tinv(4,4) = beta * (gamma - 1.0_8)
2100 
2101  ! ... compute the diagonal matrix lambda
2102  lambda(:,:) = 0.0_8
2103  do jj = 1, nd
2104  lambda(jj,jj) = ucon
2105  end do
2106  lambda(nd+1,nd+1) = ucon + spd_snd_roe * sqrt(xi_x**2 + xi_y**2)
2107  lambda(nd+2,nd+2) = ucon - spd_snd_roe * sqrt(xi_x**2 + xi_y**2)
2108 
2109  else if (nd == 3) then
2110  ! ............................................................
2111  ! ... Modified by Pooya (2018) ...............................
2112  ! ... Reference: Nonomura & Fujii, JCP (2017) ................
2113  ! ............................................................
2114  ! Ported to PlasCom2 [mtc]
2115 
2116 
2117  k_x = metricdata(1); k_y= metricdata(2); k_z = metricdata(3)
2118  ja = dsqrt(k_x**2+k_y**2+k_z**2)
2119  k_x = k_x/ja; k_y = k_y/ja; k_z = k_z/ja
2120 
2121  l_x = metricdata(4); l_y= metricdata(5); l_z = metricdata(6)
2122  ja = dsqrt(l_x**2+l_y**2+l_z**2)
2123  l_x = l_x/ja; l_y = l_y/ja; l_z = l_z/ja
2124 
2125  m_x = metricdata(7); m_y= metricdata(8); m_z = metricdata(9)
2126  ja =dsqrt(m_x**2+m_y**2+m_z**2)
2127  m_x = m_x/ja; m_y = m_y/ja; m_z = m_z/ja
2128 
2129  u2 = ux_roe**2+uy_roe**2+uz_roe**2
2130 
2131  uhat = k_x * ux_roe + k_y * uy_roe + k_z * uz_roe
2132  vhat = l_x * ux_roe + l_y * uy_roe + l_z * uz_roe
2133  what = m_x * ux_roe + m_y * uy_roe + m_z * uz_roe
2134 
2135  alpha = rho / (2.0_8 * spd_snd_roe)
2136 
2137  tmat(1,1) = 1.0_8
2138  tmat(1,2) = 0.0_8
2139  tmat(1,3) = 0.0_8
2140  tmat(1,4) = alpha
2141  tmat(1,5) = alpha
2142  tmat(2,1) = ux_roe
2143  tmat(2,2) = rho * l_x
2144  tmat(2,3) = rho * m_x
2145  tmat(2,4) = alpha * (ux_roe + spd_snd_roe * k_x)
2146  tmat(2,5) = alpha * (ux_roe - spd_snd_roe * k_x)
2147  tmat(3,1) = uy_roe
2148  tmat(3,2) = rho * l_y
2149  tmat(3,3) = rho * m_y
2150  tmat(3,4) = alpha * (uy_roe + spd_snd_roe * k_y)
2151  tmat(3,5) = alpha * (uy_roe - spd_snd_roe * k_y)
2152  tmat(4,1) = uz_roe
2153  tmat(4,2) = rho * l_z
2154  tmat(4,3) = rho * m_z
2155  tmat(4,4) = alpha * (uz_roe + spd_snd_roe * k_z)
2156  tmat(4,5) = alpha * (uz_roe - spd_snd_roe * k_z)
2157  tmat(5,1) = 0.5_8 * u2
2158  tmat(5,2) = rho * vhat
2159  tmat(5,3) = rho * what
2160  tmat(5,4) = alpha * (h_roe + spd_snd_roe * uhat)
2161  tmat(5,5) = alpha * (h_roe - spd_snd_roe * uhat)
2162 
2163  b1 = (gamma - 1.0_8) / spd_snd_roe**2
2164  b2 = 0.5_8 * (gamma - 1.0_8) * (ux_roe**2 + uy_roe**2 + uz_roe**2) / spd_snd_roe**2
2165  beta = 1.0_8 / (2.0_8 * alpha)
2166  tinv(1,1) = 1.0_8 - b2
2167  tinv(1,2) = b1 * ux_roe
2168  tinv(1,3) = b1 * uy_roe
2169  tinv(1,4) = b1 * uz_roe
2170  tinv(1,5) = -b1
2171  tinv(2,1) = -vhat/ rho
2172  tinv(2,2) = l_x/ rho
2173  tinv(2,3) = l_y/ rho
2174  tinv(2,4) = l_z/ rho
2175  tinv(2,5) = 0.0_8
2176  tinv(3,1) = -what/ rho
2177  tinv(3,2) = m_x/ rho
2178  tinv(3,3) = m_y/ rho
2179  tinv(3,4) = m_z/ rho
2180  tinv(3,5) = 0.0_8
2181  tinv(4,1) = beta * (b2 - uhat / spd_snd_roe)
2182  tinv(4,2) = -beta * (-k_x / spd_snd_roe + b1 * ux_roe)
2183  tinv(4,3) = -beta * (-k_y / spd_snd_roe + b1 * uy_roe)
2184  tinv(4,4) = -beta * (-k_z / spd_snd_roe + b1 * uz_roe)
2185  tinv(4,5) = beta * b1
2186  tinv(5,1) = beta * (b2 + uhat / spd_snd_roe)
2187  tinv(5,2) = -beta * (k_x / spd_snd_roe + b1 * ux_roe)
2188  tinv(5,3) = -beta * (k_y / spd_snd_roe + b1 * uy_roe)
2189  tinv(5,4) = -beta * (k_z / spd_snd_roe + b1 * uz_roe)
2190  tinv(5,5) = beta * b1
2191 
2192  !check=matmul(Tmat,Tinv)
2193  !do jj = 1, 5
2194  ! do i = 1, 5
2195  ! if (i==jj.and.dabs(check(i,jj)-1.0_8)>1d-11) print*,dir,i,jj,check(i,jj)
2196  ! if (i.ne.jj.and.dabs(check(i,jj))>1d-11) print*,dir,i,jj,check(i,jj)
2197  ! enddo
2198  !enddo
2199 
2200  ! ... compute the diagonal matrix lambda
2201  lambda(:,:) = 0.0_8
2202  do jj = 1, nd
2203  lambda(jj,jj) = ucon
2204  end do
2205  lambda(nd+1,nd+1) = ucon + spd_snd_roe * dsqrt(xi_x**2 + xi_y**2 + xi_z**2)
2206  lambda(nd+2,nd+2) = ucon - spd_snd_roe * dsqrt(xi_x**2 + xi_y**2 + xi_z**2)
2207 
2208  end if
2209 
2210  end subroutine sat_form_roe_matrices_test
2211 
2212  SUBROUTINE dissipationweight(numDim,bufferSize,numPoints,bufferInterval, &
2213  sigmaDissipation,sigmaDilatation,dilatationRamp,dilatationCutoff, &
2214  divV,sigmaDiss)
2215 
2216  IMPLICIT NONE
2217 
2218  INTEGER(KIND=4), INTENT(IN) :: numDim
2219  INTEGER(KIND=8), INTENT(IN) :: bufferSize(numdim), numPoints
2220  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
2221  REAL(KIND=8), INTENT(IN) :: sigmaDissipation,sigmaDilatation
2222  REAL(KIND=8), INTENT(IN) :: dilatationRamp,dilatationCutoff
2223  REAL(KIND=8), INTENT(IN) :: divV(numpoints)
2224  REAL(KIND=8), INTENT(OUT) :: sigmaDiss(numpoints)
2225 
2226  INTEGER(KIND=8) :: I, J, K
2227  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
2228  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
2229 
2230  istart = bufferinterval(1)
2231  iend = bufferinterval(2)
2232  xsize = buffersize(1)
2233 
2234  IF(numdim == 1) THEN
2235  DO i = istart, iend
2236  sigmadiss(i) = (sigmadissipation + &
2237  sigmadilatation*0.5_8*(1.0_8 + tanh(dilatationramp*(dilatationcutoff-divv(i)))))
2238  END DO
2239  ELSE IF(numdim == 2) THEN
2240  jstart = bufferinterval(3)
2241  jend = bufferinterval(4)
2242  DO j = jstart, jend
2243  yindex = (j-1)*xsize
2244  DO i = istart, iend
2245  bufferindex = yindex + i
2246  sigmadiss(bufferindex) = (sigmadissipation + &
2247  sigmadilatation*0.5_8*(1.0_8 + tanh(dilatationramp*(dilatationcutoff-divv(bufferindex)))))
2248  END DO
2249  END DO
2250  ELSE IF(numdim == 3) THEN
2251  jstart = bufferinterval(3)
2252  jend = bufferinterval(4)
2253  kstart = bufferinterval(5)
2254  kend = bufferinterval(6)
2255  nplane = xsize * buffersize(2)
2256  DO k = kstart, kend
2257  zindex = (k-1)*nplane
2258  DO j = jstart, jend
2259  yzindex = zindex + (j-1)*xsize
2260  DO i = istart, iend
2261  bufferindex = yzindex + i
2262  sigmadiss(bufferindex) = (sigmadissipation + &
2263  sigmadilatation*0.5_8*(1.0_8 + tanh(dilatationramp*(dilatationcutoff-divv(bufferindex)))))
2264  END DO
2265  END DO
2266  END DO
2267  ENDIF
2268 
2269 
2270  END SUBROUTINE dissipationweight
2271 
2272  subroutine sat_form_roe_matrices_modified(u_in, u_out, gamma, norm, MT1, dir, ND, tmat, tinv, lambda)
2273 
2274 ! USE ModGlobal
2275 
2276  Implicit None
2277 
2278  Real(8), Pointer :: u_in(:), u_out(:), tmat(:,:), tinv(:,:), norm(:), lambda(:,:), MT1(:)
2279  Real(8) :: gamma, rho_in, ux_in, uy_in, uz_in, ke_in, p_in, en_in, h_in
2280  Real(8) :: gm1, rho_out, ux_out, uy_out, uz_out, ke_out, p_out, en_out, h_out
2281  Real(8) :: cc, bb, ux_Roe, uy_Roe, uz_Roe, h_Roe, ke_Roe, spd_snd_Roe, rho_Roe
2282  Real(8) :: XI_X_TILDE, XI_Y_TILDE, XI_Z_TILDE, ucon, rho, alpha, theta, beta, phi2
2283  Real(8) :: XI_X, XI_Y, XI_Z, denom, b1, b2, U2
2284  Integer :: N, ND, jj, dir, i
2285  Real(8) :: uhat, vhat, what, k_x, k_y, k_z, l_x, l_y, l_z, m_x, m_y, m_z, ja!,check(5,5)
2286 
2287 
2288 !NEW
2289  ! ... constants
2290  gm1 = gamma - 1.0_8
2291 
2292  ! ... form primative variables from "in"
2293  rho_in = u_in(1)
2294  ux_in = u_in(2) / rho_in
2295  uy_in = 0.0_8
2296  uz_in = 0.0_8
2297  if (nd >= 2) uy_in = u_in(3) / rho_in
2298  if (nd == 3) uz_in = u_in(4) / rho_in
2299 
2300  ke_in = 0.5_8 * (ux_in**2 + uy_in**2 + uz_in**2)
2301  en_in = u_in(nd+2)/rho_in - ke_in
2302  p_in = gm1 * rho_in * en_in
2303  h_in = (gamma / gm1) * p_in / rho_in + ke_in
2304  ! ... form primative variables from "out"
2305  rho_out = u_out(1)
2306  ux_out = u_out(2) / rho_out
2307  uy_out = 0.0_8
2308  uz_out = 0.0_8
2309  if (nd >= 2) uy_out = u_out(3) / rho_out
2310  if (nd == 3) uz_out = u_out(4) / rho_out
2311 
2312  ke_out = 0.5_8 * (ux_out**2 + uy_out**2 + uz_out**2)
2313  en_out = u_out(nd+2)/rho_out - ke_out
2314  p_out = gm1 * rho_out * en_out
2315  h_out = (gamma / gm1) * p_out / rho_out + ke_out
2316 
2317  ! ... form Roe variables
2318  cc = sqrt(rho_out / rho_in)
2319  bb = 1.0_8 / (1.0_8 + cc)
2320  ux_roe = (ux_in + ux_out * cc) * bb
2321  uy_roe = (uy_in + uy_out * cc) * bb
2322  uz_roe = (uz_in + uz_out * cc) * bb
2323  h_roe = ( h_in + h_out * cc) * bb
2324  ke_roe = 0.5_8 * (ux_roe**2 + uy_roe**2 + uz_roe**2)
2325 
2326  if (h_roe - ke_roe <= 0.0_8) then
2327  write (*,'(5(a,E20.8,1X))') "gamma = ", gamma, ", ke_in = ", ke_in, ", en_in = ", en_in, &
2328  ", p_in = ", p_in, ", h_in = ", h_in
2329  stop
2330  end if
2331 
2332  spd_snd_roe = sqrt(gm1*(h_roe - ke_roe))
2333  rho_roe = rho_in * cc
2334 
2335  ! ... variable renaming
2336  xi_x = norm(1); xi_y = 0.0_8; xi_z = 0.0_8
2337  if (nd >= 2) xi_y = norm(2)
2338  if (nd == 3) xi_z = norm(3)
2339  denom = 1.0_8 / sqrt(xi_x**2 + xi_y**2 + xi_z**2)
2340  xi_x_tilde = xi_x * denom
2341  xi_y_tilde = xi_y * denom
2342  xi_z_tilde = xi_z * denom
2343  rho = rho_roe
2344 
2345  ! ... form transformation matrices
2346  ucon = ux_roe * xi_x + uy_roe * xi_y + uz_roe * xi_z
2347 
2348 
2349  if (nd == 2) then
2350 
2351  ! ... directly from Pulliam & Chaussee
2352  alpha = 0.5_8 * rho / spd_snd_roe
2353  beta = 1.0_8 / (rho * spd_snd_roe)
2354  alpha = rho / (sqrt(2.0_8) * spd_snd_roe)
2355  beta = 1.0_8 / (sqrt(2.0_8) * rho * spd_snd_roe)
2356  theta = xi_x_tilde * ux_roe + xi_y_tilde * uy_roe
2357  phi2 = 0.5_8 * (gamma - 1.0_8) * (ux_roe**2 + uy_roe**2)
2358 
2359  tmat(1,1) = 1.0_8
2360  tmat(1,2) = 0.0_8
2361  tmat(1,3) = alpha
2362  tmat(1,4) = alpha
2363  tmat(2,1) = ux_roe
2364  tmat(2,2) = xi_y_tilde * rho
2365  tmat(2,3) = alpha * (ux_roe + xi_x_tilde * spd_snd_roe)
2366  tmat(2,4) = alpha * (ux_roe - xi_x_tilde * spd_snd_roe)
2367  tmat(3,1) = uy_roe
2368  tmat(3,2) = -xi_x_tilde * rho
2369  tmat(3,3) = alpha * (uy_roe + xi_y_tilde * spd_snd_roe)
2370  tmat(3,4) = alpha * (uy_roe - xi_y_tilde * spd_snd_roe)
2371  tmat(4,1) = phi2 / (gamma - 1.0_8)
2372  tmat(4,2) = rho * (xi_y_tilde * ux_roe - xi_x_tilde * uy_roe)
2373  tmat(4,3) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) + spd_snd_roe * theta)
2374  tmat(4,4) = alpha * ((phi2 + spd_snd_roe**2)/(gamma - 1.0_8) - spd_snd_roe * theta)
2375 
2376  tinv(1,1) = 1.0_8 - phi2 / spd_snd_roe**2
2377  tinv(1,2) = (gamma - 1.0_8) * ux_roe / spd_snd_roe**2
2378  tinv(1,3) = (gamma - 1.0_8) * uy_roe / spd_snd_roe**2
2379  tinv(1,4) = -(gamma - 1.0_8) / spd_snd_roe**2
2380  tinv(2,1) = -(xi_y_tilde * ux_roe - xi_x_tilde * uy_roe) / rho
2381  tinv(2,2) = xi_y_tilde / rho
2382  tinv(2,3) = -xi_x_tilde / rho
2383  tinv(2,4) = 0.0_8
2384  tinv(3,1) = beta * (phi2 - spd_snd_roe * theta)
2385  tinv(3,2) = beta * (xi_x_tilde * spd_snd_roe - (gamma - 1.0_8) * ux_roe)
2386  tinv(3,3) = beta * (xi_y_tilde * spd_snd_roe - (gamma - 1.0_8) * uy_roe)
2387  tinv(3,4) = beta * (gamma - 1.0_8)
2388  tinv(4,1) = beta * (phi2 + spd_snd_roe * theta)
2389  tinv(4,2) = -beta * (xi_x_tilde * spd_snd_roe + (gamma - 1.0_8) * ux_roe)
2390  tinv(4,3) = -beta * (xi_y_tilde * spd_snd_roe + (gamma - 1.0_8) * uy_roe)
2391  tinv(4,4) = beta * (gamma - 1.0_8)
2392 
2393  ! ... compute the diagonal matrix lambda
2394  lambda(:,:) = 0.0_8
2395  do jj = 1, nd
2396  lambda(jj,jj) = ucon
2397  end do
2398  lambda(nd+1,nd+1) = ucon + spd_snd_roe * sqrt(xi_x**2 + xi_y**2)
2399  lambda(nd+2,nd+2) = ucon - spd_snd_roe * sqrt(xi_x**2 + xi_y**2)
2400 ! --- NEW ---
2401 
2402  else if (nd == 3) then
2403  ! ............................................................
2404  ! ... Modified by Pooya (2018) ...............................
2405  ! ... Reference: Nonomura & Fujii, JCP (2017) ................
2406  ! ............................................................
2407  ! Ported to PlasCom2 [mtc]
2408 
2409 ! if (dir==1) then
2410 
2411  k_x = mt1(1); k_y= mt1(2); k_z = mt1(3)
2412  ja = dsqrt(k_x**2+k_y**2+k_z**2)
2413  k_x = k_x/ja; k_y = k_y/ja; k_z = k_z/ja
2414 
2415  l_x = mt1(4); l_y= mt1(5); l_z = mt1(6)
2416  ja = dsqrt(l_x**2+l_y**2+l_z**2)
2417  l_x = l_x/ja; l_y = l_y/ja; l_z = l_z/ja
2418 
2419  m_x = mt1(7); m_y= mt1(8); m_z = mt1(9)
2420  ja =dsqrt(m_x**2+m_y**2+m_z**2)
2421  m_x = m_x/ja; m_y = m_y/ja; m_z = m_z/ja
2422 
2423 ! elseif (dir==2) then
2424 !
2425 ! l_x = MT1(1); l_y= MT1(4); l_z = MT1(7)
2426 ! ja = dsqrt(l_x**2+l_y**2+l_z**2)
2427 ! l_x = l_x/ja; l_y = l_y/ja; l_z = l_z/ja
2428 !
2429 ! k_x = MT1(2); k_y= MT1(5); k_z = MT1(8)
2430 ! ja = dsqrt(k_x**2+k_y**2+k_z**2)
2431 ! k_x = k_x/ja; k_y = k_y/ja; k_z = k_z/ja
2432 !
2433 ! m_x = MT1(3); m_y= MT1(6); m_z = MT1(9)
2434 ! ja =dsqrt(m_x**2+m_y**2+m_z**2)
2435 ! m_x = m_x/ja; m_y = m_y/ja; m_z = m_z/ja
2436 !
2437 ! elseif (dir==3) then
2438 !
2439 ! m_x = MT1(1); m_y= MT1(4); m_z = MT1(7)
2440 ! ja = dsqrt(m_x**2+m_y**2+m_z**2)
2441 ! m_x = m_x/ja; m_y = m_y/ja; m_z = m_z/ja
2442 !
2443 ! l_x = MT1(2); l_y= MT1(5); l_z = MT1(8)
2444 ! ja = dsqrt(l_x**2+l_y**2+l_z**2)
2445 ! l_x = l_x/ja; l_y = l_y/ja; l_z = l_z/ja
2446 !
2447 ! k_x = MT1(3); k_y= MT1(6); k_z = MT1(9)
2448 ! ja =dsqrt(k_x**2+k_y**2+k_z**2)
2449 ! k_x = k_x/ja; k_y = k_y/ja; k_z = k_z/ja
2450 !
2451 !
2452 
2453  u2 = ux_roe**2+uy_roe**2+uz_roe**2
2454 
2455  uhat = k_x * ux_roe + k_y * uy_roe + k_z * uz_roe
2456  vhat = l_x * ux_roe + l_y * uy_roe + l_z * uz_roe
2457  what = m_x * ux_roe + m_y * uy_roe + m_z * uz_roe
2458 
2459  alpha = rho / (2.0_8 * spd_snd_roe)
2460 
2461  tmat(1,1) = 1.0_8
2462  tmat(1,2) = 0.0_8
2463  tmat(1,3) = 0.0_8
2464  tmat(1,4) = alpha
2465  tmat(1,5) = alpha
2466  tmat(2,1) = ux_roe
2467  tmat(2,2) = rho * l_x
2468  tmat(2,3) = rho * m_x
2469  tmat(2,4) = alpha * (ux_roe + spd_snd_roe * k_x)
2470  tmat(2,5) = alpha * (ux_roe - spd_snd_roe * k_x)
2471  tmat(3,1) = uy_roe
2472  tmat(3,2) = rho * l_y
2473  tmat(3,3) = rho * m_y
2474  tmat(3,4) = alpha * (uy_roe + spd_snd_roe * k_y)
2475  tmat(3,5) = alpha * (uy_roe - spd_snd_roe * k_y)
2476  tmat(4,1) = uz_roe
2477  tmat(4,2) = rho * l_z
2478  tmat(4,3) = rho * m_z
2479  tmat(4,4) = alpha * (uz_roe + spd_snd_roe * k_z)
2480  tmat(4,5) = alpha * (uz_roe - spd_snd_roe * k_z)
2481  tmat(5,1) = 0.5_8 * u2
2482  tmat(5,2) = rho * vhat
2483  tmat(5,3) = rho * what
2484  tmat(5,4) = alpha * (h_roe + spd_snd_roe * uhat)
2485  tmat(5,5) = alpha * (h_roe - spd_snd_roe * uhat)
2486 
2487  b1 = (gamma - 1.0_8) / spd_snd_roe**2
2488  b2 = 0.5_8 * (gamma - 1.0_8) * (ux_roe**2 + uy_roe**2 + uz_roe**2) / spd_snd_roe**2
2489  beta = 1.0_8 / (2.0_8 * alpha)
2490  tinv(1,1) = 1.0_8 - b2
2491  tinv(1,2) = b1 * ux_roe
2492  tinv(1,3) = b1 * uy_roe
2493  tinv(1,4) = b1 * uz_roe
2494  tinv(1,5) = -b1
2495  tinv(2,1) = -vhat/ rho
2496  tinv(2,2) = l_x/ rho
2497  tinv(2,3) = l_y/ rho
2498  tinv(2,4) = l_z/ rho
2499  tinv(2,5) = 0.0_8
2500  tinv(3,1) = -what/ rho
2501  tinv(3,2) = m_x/ rho
2502  tinv(3,3) = m_y/ rho
2503  tinv(3,4) = m_z/ rho
2504  tinv(3,5) = 0.0_8
2505  tinv(4,1) = beta * (b2 - uhat / spd_snd_roe)
2506  tinv(4,2) = -beta * (-k_x / spd_snd_roe + b1 * ux_roe)
2507  tinv(4,3) = -beta * (-k_y / spd_snd_roe + b1 * uy_roe)
2508  tinv(4,4) = -beta * (-k_z / spd_snd_roe + b1 * uz_roe)
2509  tinv(4,5) = beta * b1
2510  tinv(5,1) = beta * (b2 + uhat / spd_snd_roe)
2511  tinv(5,2) = -beta * (k_x / spd_snd_roe + b1 * ux_roe)
2512  tinv(5,3) = -beta * (k_y / spd_snd_roe + b1 * uy_roe)
2513  tinv(5,4) = -beta * (k_z / spd_snd_roe + b1 * uz_roe)
2514  tinv(5,5) = beta * b1
2515 
2516  !check=matmul(Tmat,Tinv)
2517  !do jj = 1, 5
2518  ! do i = 1, 5
2519  ! if (i==jj.and.dabs(check(i,jj)-1.0_8)>1d-11) print*,dir,i,jj,check(i,jj)
2520  ! if (i.ne.jj.and.dabs(check(i,jj))>1d-11) print*,dir,i,jj,check(i,jj)
2521  ! enddo
2522  !enddo
2523 
2524  ! ... compute the diagonal matrix lambda
2525  lambda(:,:) = 0.0_8
2526  do jj = 1, nd
2527  lambda(jj,jj) = ucon
2528  end do
2529  lambda(nd+1,nd+1) = ucon + spd_snd_roe * dsqrt(xi_x**2 + xi_y**2 + xi_z**2)
2530  lambda(nd+2,nd+2) = ucon - spd_snd_roe * dsqrt(xi_x**2 + xi_y**2 + xi_z**2)
2531 
2532  end if
2533 
2534  end subroutine sat_form_roe_matrices_modified
2535 
2536  subroutine sat_form_roe_matrices_multispecies(u_in, u_out, p_in, p_out, MT1, JAC, u_ref, ND, dir, &
2537  nAuxVars, gamma, tmat, tinv, lambda)
2538 
2539 ! USE ModGlobal
2540 ! USE ModDataStruct
2541 
2542  Implicit None
2543 
2544  Real(8), Pointer :: u_in(:), u_out(:), tmat(:,:), tinv(:,:), lambda(:,:), MT1(:)
2545  Real(8) :: gamma, rho_in, ux_in, uy_in, uz_in, ke_in, p_in, en_in, H_in, JAC
2546  Real(8) :: gm1, rho_out, ux_out, uy_out, uz_out, ke_out, p_out, en_out, H_out
2547  Real(8) :: cc, bb, ux_Roe, uy_Roe, uz_Roe, h_Roe, en_Roe, ke_Roe, spd_snd_Roe, rho_Roe
2548  Real(8) :: XI_X_TILDE, XI_Y_TILDE, XI_Z_TILDE, rho, alpha, theta, beta, phi2
2549  Real(8) :: XI_X, XI_Y, XI_Z, denom
2550  Integer :: N,ND,jj,nAuxVars,dir
2551 ! type(t_mixt), pointer :: state
2552  Real(8) :: Y_in(nauxvars),Y_out(nauxvars),ja,Residual,dp_drho_ROE,dp_de_ROE,dp_drhoY_ROE(nauxvars)
2553  Real(8) :: z_Roe(nauxvars),Y_ROE(nauxvars),u_ref,b1,b2,b3,dp_drho_bar,dp_de_bar,dp_drhoY_bar(nauxvars)
2554  Real(8) :: metric_vec(3),p_ROE,temp_ROE,junk
2555  real(8), pointer :: pnorm_vec(:)
2556  Real(8) :: uhat, vhat, what, k_x, k_y, k_z, l_x, l_y, l_z, m_x, m_y, m_z!, check(8,8)
2557 ! integer :: i
2558 
2559  If (nauxvars == 0) then
2560 
2561  nullify(pnorm_vec)
2562  allocate(pnorm_vec(nd))
2563  xi_x = mt1(dir) * jac
2564  xi_y = mt1(dir+nd) * jac
2565  xi_z = 0.0_8
2566  if (nd >= 3) then
2567  xi_z = mt1(dir+2*nd) *jac
2568  end if
2569  metric_vec = (/ xi_x, xi_y, xi_z /)
2570  pnorm_vec(:) = metric_vec(1:nd)
2571 ! Call SAT_Form_Roe_Matrices(u_in, u_out, gamma, pnorm_vec, Tmat, Tinv, Lambda)
2573  Call sat_form_roe_matrices_modified(u_in, u_out, gamma, pnorm_vec, mt1, dir, nd, tmat, tinv, lambda)
2574  deallocate(pnorm_vec)
2575  nullify(pnorm_vec)
2576 
2577  elseif (nauxvars > 0) then
2578  ! ... form primative variables from "in"
2579  rho_in = u_in(1)
2580  ux_in = u_in(2) / rho_in
2581  uy_in = u_in(3) / rho_in
2582  uz_in = 0.0_8
2583  if (nd == 3) uz_in = u_in(4) / rho_in
2584  ke_in = 0.5_8 * (ux_in**2 + uy_in**2 + uz_in**2)
2585  en_in = u_in(nd+2)/rho_in - ke_in
2586  h_in = (u_in(nd+2)+p_in) / rho_in
2587  p_in = p_in
2588  do jj = 1, nauxvars
2589  y_in(jj) = u_in(nd+2+jj)/ rho_in
2590  enddo
2591 
2592  ! ... form primative variables from "out"
2593  rho_out = u_out(1)
2594  ux_out = u_out(2) / rho_out
2595  uy_out = u_out(3) / rho_out
2596  uz_out = 0.0_8
2597  if (nd == 3) uz_out = u_out(4) / rho_out
2598  ke_out = 0.5_8 * (ux_out**2 + uy_out**2 + uz_out**2)
2599  en_out = u_out(nd+2)/rho_out - ke_out
2600  h_out = (u_out(nd+2)+p_out) / rho_out
2601  p_out = p_out
2602  do jj = 1, nauxvars
2603  y_out(jj) = u_out(nd+2+jj)/ rho_out
2604  enddo
2605 
2606  ! ... form Roe variables
2607  cc = sqrt(rho_out / rho_in)
2608  bb = 1.0_8 / (1.0_8 + cc)
2609  ux_roe = (ux_in + ux_out * cc) * bb
2610  uy_roe = (uy_in + uy_out * cc) * bb
2611  uz_roe = (uz_in + uz_out * cc) * bb
2612  h_roe = ( h_in + h_out * cc) * bb
2613  en_roe = ( en_in + en_out * cc) * bb
2614  ke_roe = 0.5_8 * (ux_roe**2 + uy_roe**2 + uz_roe**2)
2615  rho_roe = rho_in * cc
2616  do jj = 1, nauxvars
2617  y_roe(jj) = (y_in(jj) + y_out(jj) * cc) * bb
2618  enddo
2619 
2620  ! ... variable renaming
2621  xi_x = mt1(dir) * jac
2622  xi_y = mt1(dir+nd) * jac
2623  xi_z = 0.0_8
2624  if (nd >= 3) then
2625  xi_z = mt1(dir+2*nd) *jac
2626  end if
2627  denom = 1.0_8 / dsqrt(xi_x**2 + xi_y**2 + xi_z**2)
2628  xi_x_tilde = xi_x * denom
2629  xi_y_tilde = xi_y * denom
2630  xi_z_tilde = xi_z * denom
2631  rho = rho_roe
2632 
2633 ! CALL PDerivative_ROE(RHO,en_ROE,Y_ROE,u_ref,nAuxvars,state,dp_drho_bar,dp_de_bar,dp_drhoY_bar)
2634  CALL pderivative_roe(rho,en_roe,y_roe,u_ref,nauxvars,dp_drho_bar,dp_de_bar,dp_drhoy_bar)
2635 
2636  residual = (p_out-p_in)-(dp_drho_bar*(rho_out-rho_in)+dp_de_bar*(en_out-en_in))
2637  do jj = 1, nauxvars
2638  residual = residual - dp_drhoy_bar(jj) * (u_out(nd+2+jj)-u_in(nd+2+jj))
2639  enddo
2640 
2641  ja = (dp_de_bar*(en_out-en_in))**2 + (dp_drho_bar*(rho_out-rho_in))**2
2642  do jj = 1, nauxvars
2643  ja = ja + (dp_drhoy_bar(jj)*(u_out(nd+2+jj)-u_in(nd+2+jj)))**2
2644  enddo
2645 
2646  If (residual == 0.0_8.or.ja<1d-13) then
2647  dp_de_roe = dp_de_bar
2648  dp_drho_roe = dp_drho_bar
2649  dp_drhoy_roe = dp_drhoy_bar
2650  else
2651  dp_de_roe = dp_de_bar * (1.0_8+dp_de_bar*(en_out-en_in)/ja*residual)
2652  dp_drho_roe = dp_drho_bar * (1.0_8+dp_drho_bar*(rho_out-rho_in)/ja*residual)
2653  do jj =1 , nauxvars
2654  dp_drhoy_roe(jj) = dp_drhoy_bar(jj) * (1.0_8+dp_drhoy_bar(jj)*(u_out(nd+2+jj)-u_in(nd+2+jj))/ja*residual)
2655  enddo
2656  endif
2657 
2658  do jj =1, nauxvars
2659  z_roe(jj) = -rho * dp_drhoy_roe(jj)/ dp_de_roe
2660  enddo
2661 
2662  spd_snd_roe = dp_drho_roe + dp_de_roe/ rho * (h_roe -en_roe -ke_roe)
2663  do jj = 1, nauxvars
2664  spd_snd_roe = spd_snd_roe + y_roe(jj) * dp_drhoy_roe(jj)
2665  enddo
2666  spd_snd_roe = dsqrt(spd_snd_roe)
2667 
2668  b1 = dp_de_roe / rho / spd_snd_roe**2
2669  b2 = 1.0_8 + b1 * (ux_roe**2 + uy_roe**2 + uz_roe**2) - b1 * h_roe
2670  b3 = 0.0_8
2671  do jj = 1, nauxvars
2672  b3 = b3 + y_roe(jj)* z_roe(jj)
2673  enddo
2674  b3 = b3 * b1
2675 
2676  theta = xi_x_tilde * ux_roe + xi_y_tilde * uy_roe + xi_z_tilde * uz_roe
2677  tmat =0.0_8; tinv = 0.0_8
2678 
2679  if (nd == 2) then
2680 
2681  alpha = rho / (sqrt(2.0_8) * spd_snd_roe)
2682  beta = 1.0_8 / (sqrt(2.0_8) * rho * spd_snd_roe)
2683 
2684  tmat(1,1) = 1.0_8
2685  tmat(1,2) = 0.0_8
2686  tmat(1,3) = alpha
2687  tmat(1,4) = alpha
2688  tmat(2,1) = ux_roe
2689  tmat(2,2) = xi_y_tilde * rho
2690  tmat(2,3) = alpha * (ux_roe + xi_x_tilde * spd_snd_roe)
2691  tmat(2,4) = alpha * (ux_roe - xi_x_tilde * spd_snd_roe)
2692  tmat(3,1) = uy_roe
2693  tmat(3,2) = -xi_x_tilde * rho
2694  tmat(3,3) = alpha * (uy_roe + xi_y_tilde * spd_snd_roe)
2695  tmat(3,4) = alpha * (uy_roe - xi_y_tilde * spd_snd_roe)
2696  tmat(4,1) = h_roe - rho * spd_snd_roe**2 / dp_de_roe
2697  tmat(4,2) = rho * (xi_y_tilde * ux_roe - xi_x_tilde * uy_roe)
2698  tmat(4,3) = alpha * (h_roe + spd_snd_roe * theta)
2699  tmat(4,4) = alpha * (h_roe - spd_snd_roe * theta)
2700 
2701  do jj = 1, nauxvars
2702  tmat(4,4+jj) = z_roe(jj)
2703  tmat(4+jj,1) = y_roe(jj)
2704  tmat(4+jj,3) = alpha * y_roe(jj)
2705  tmat(4+jj,4) = alpha * y_roe(jj)
2706  tmat(4+jj,4+jj) = 1.0_8
2707  enddo
2708 
2709  tinv(1,1) = 1.0_8 - b2 - b3
2710  tinv(1,2) = b1 * ux_roe
2711  tinv(1,3) = b1 * uy_roe
2712  tinv(1,4) = -b1
2713  tinv(2,1) = -(xi_y_tilde * ux_roe - xi_x_tilde * uy_roe) / rho
2714  tinv(2,2) = xi_y_tilde / rho
2715  tinv(2,3) = -xi_x_tilde / rho
2716  tinv(3,1) = beta * spd_snd_roe**2 * (b2 + b3 - theta / spd_snd_roe)
2717  tinv(3,2) = beta * spd_snd_roe * (xi_x_tilde - b1 * ux_roe * spd_snd_roe)
2718  tinv(3,3) = beta * spd_snd_roe * (xi_y_tilde - b1 * uy_roe * spd_snd_roe)
2719  tinv(3,4) = beta * b1 * spd_snd_roe**2
2720  tinv(4,1) = beta * spd_snd_roe**2 * (b2 + b3 + theta / spd_snd_roe)
2721  tinv(4,2) = -beta * spd_snd_roe * (xi_x_tilde + b1 * ux_roe * spd_snd_roe)
2722  tinv(4,3) = -beta * spd_snd_roe * (xi_y_tilde + b1 * uy_roe * spd_snd_roe)
2723  tinv(4,4) = beta * b1 * spd_snd_roe**2
2724 
2725  do jj = 1, nauxvars
2726  tinv(1,4+jj) = b1 * z_roe(jj)
2727  tinv(3,4+jj) = -beta * b1 * spd_snd_roe**2 * z_roe(jj)
2728  tinv(4,4+jj) = -beta * b1 * spd_snd_roe**2 * z_roe(jj)
2729  tinv(4+jj,1) = -y_roe(jj)
2730  tinv(4+jj,4+jj) = 1.0_8
2731  enddo
2732 
2733  ! ... compute the diagonal matrix lambda
2734  lambda(:,:) = 0.0_8
2735  do jj = 1, nd
2736  lambda(jj,jj) = theta
2737  end do
2738  lambda(nd+1,nd+1) = theta + spd_snd_roe * dsqrt(xi_x**2 + xi_y**2)
2739  lambda(nd+2,nd+2) = theta - spd_snd_roe * dsqrt(xi_x**2 + xi_y**2)
2740  do jj = 1, nauxvars
2741  lambda(4+jj,4+jj) = theta
2742  enddo
2743 
2744  elseif (nd ==3) then
2745  ! ............................................................
2746  ! ... Modified by Pooya (2018) ...............................
2747  ! ... Reference: Nonomura & Fujii, JCP (2017) ................
2748  ! ............................................................
2749  if (dir==1) then
2750  k_x = mt1(1)*jac; k_y= mt1(4)*jac; k_z = mt1(7)*jac
2751  ja = dsqrt(k_x**2+k_y**2+k_z**2)
2752  k_x = k_x/ja; k_y = k_y/ja; k_z = k_z/ja
2753 
2754  l_x = mt1(2)*jac; l_y= mt1(5)*jac; l_z = mt1(8)*jac
2755  ja = dsqrt(l_x**2+l_y**2+l_z**2)
2756  l_x = l_x/ja; l_y = l_y/ja; l_z = l_z/ja
2757 
2758  m_x = mt1(3)*jac; m_y= mt1(6)*jac; m_z = mt1(9)*jac
2759  ja =dsqrt(m_x**2+m_y**2+m_z**2)
2760  m_x = m_x/ja; m_y = m_y/ja; m_z = m_z/ja
2761 
2762  elseif (dir==2) then
2763  l_x = mt1(1)*jac; l_y= mt1(4)*jac; l_z = mt1(7)*jac
2764  ja = dsqrt(l_x**2+l_y**2+l_z**2)
2765  l_x = l_x/ja; l_y = l_y/ja; l_z = l_z/ja
2766 
2767  k_x = mt1(2)*jac; k_y= mt1(5)*jac; k_z = mt1(8)*jac
2768  ja = dsqrt(k_x**2+k_y**2+k_z**2)
2769  k_x = k_x/ja; k_y = k_y/ja; k_z = k_z/ja
2770 
2771  m_x = mt1(3)*jac; m_y= mt1(6)*jac; m_z = mt1(9)*jac
2772  ja =dsqrt(m_x**2+m_y**2+m_z**2)
2773  m_x = m_x/ja; m_y = m_y/ja; m_z = m_z/ja
2774 
2775  elseif (dir==3) then
2776  m_x = mt1(1)*jac; m_y= mt1(4)*jac; m_z = mt1(7)*jac
2777  ja = dsqrt(m_x**2+m_y**2+m_z**2)
2778  m_x = m_x/ja; m_y = m_y/ja; m_z = m_z/ja
2779 
2780  l_x = mt1(2)*jac; l_y= mt1(5)*jac; l_z = mt1(8)*jac
2781  ja = dsqrt(l_x**2+l_y**2+l_z**2)
2782  l_x = l_x/ja; l_y = l_y/ja; l_z = l_z/ja
2783 
2784  k_x = mt1(3)*jac; k_y= mt1(6)*jac; k_z = mt1(9)*jac
2785  ja =dsqrt(k_x**2+k_y**2+k_z**2)
2786  k_x = k_x/ja; k_y = k_y/ja; k_z = k_z/ja
2787 
2788  endif
2789 
2790  uhat = k_x * ux_roe + k_y * uy_roe + k_z * uz_roe
2791  vhat = l_x * ux_roe + l_y * uy_roe + l_z * uz_roe
2792  what = m_x * ux_roe + m_y * uy_roe + m_z * uz_roe
2793 
2794  alpha = rho / (2.0_8 * spd_snd_roe)
2795 
2796  tmat(1,1) = 1.0_8
2797  tmat(1,2) = 0.0_8
2798  tmat(1,3) = 0.0_8
2799  tmat(1,4) = alpha
2800  tmat(1,5) = alpha
2801  tmat(2,1) = ux_roe
2802  tmat(2,2) = rho * l_x
2803  tmat(2,3) = rho * m_x
2804  tmat(2,4) = alpha * (ux_roe + spd_snd_roe * k_x)
2805  tmat(2,5) = alpha * (ux_roe - spd_snd_roe * k_x)
2806  tmat(3,1) = uy_roe
2807  tmat(3,2) = rho * l_y
2808  tmat(3,3) = rho * m_y
2809  tmat(3,4) = alpha * (uy_roe + spd_snd_roe * k_y)
2810  tmat(3,5) = alpha * (uy_roe - spd_snd_roe * k_y)
2811  tmat(4,1) = uz_roe
2812  tmat(4,2) = rho * l_z
2813  tmat(4,3) = rho * m_z
2814  tmat(4,4) = alpha * (uz_roe + spd_snd_roe * k_z)
2815  tmat(4,5) = alpha * (uz_roe - spd_snd_roe * k_z)
2816  tmat(5,1) = h_roe - 1.0_8 / b1
2817  tmat(5,2) = rho * vhat
2818  tmat(5,3) = rho * what
2819  tmat(5,4) = alpha * (h_roe + spd_snd_roe * uhat)
2820  tmat(5,5) = alpha * (h_roe - spd_snd_roe * uhat)
2821 
2822  do jj = 1, nauxvars
2823  tmat(5,5+jj) = z_roe(jj)
2824  tmat(5+jj,1) = y_roe(jj)
2825  tmat(5+jj,4) = alpha * y_roe(jj)
2826  tmat(5+jj,5) = alpha * y_roe(jj)
2827  tmat(5+jj,5+jj) = 1.0_8
2828  enddo
2829 
2830  beta = 1.0_8 / (2.0_8 * alpha)
2831  tinv(1,1) = 1.0_8 - b2 -b3
2832  tinv(1,2) = b1 * ux_roe
2833  tinv(1,3) = b1 * uy_roe
2834  tinv(1,4) = b1 * uz_roe
2835  tinv(1,5) = -b1
2836  tinv(2,1) = -vhat/ rho
2837  tinv(2,2) = l_x/ rho
2838  tinv(2,3) = l_y/ rho
2839  tinv(2,4) = l_z/ rho
2840  tinv(2,5) = 0.0_8
2841  tinv(3,1) = -what/ rho
2842  tinv(3,2) = m_x/ rho
2843  tinv(3,3) = m_y/ rho
2844  tinv(3,4) = m_z/ rho
2845  tinv(3,5) = 0.0_8
2846  tinv(4,1) = beta * (b2 + b3 - uhat / spd_snd_roe)
2847  tinv(4,2) = -beta * (-k_x / spd_snd_roe + b1 * ux_roe)
2848  tinv(4,3) = -beta * (-k_y / spd_snd_roe + b1 * uy_roe)
2849  tinv(4,4) = -beta * (-k_z / spd_snd_roe + b1 * uz_roe)
2850  tinv(4,5) = beta * b1
2851  tinv(5,1) = beta * (b2 + b3 + uhat / spd_snd_roe)
2852  tinv(5,2) = -beta * (k_x / spd_snd_roe + b1 * ux_roe)
2853  tinv(5,3) = -beta * (k_y / spd_snd_roe + b1 * uy_roe)
2854  tinv(5,4) = -beta * (k_z / spd_snd_roe + b1 * uz_roe)
2855  tinv(5,5) = beta * b1
2856 
2857  do jj = 1, nauxvars
2858  tinv(1,5+jj) = b1 * z_roe(jj)
2859  tinv(4,5+jj) = -b1 * z_roe(jj) * spd_snd_roe / rho
2860  tinv(5,5+jj) = -b1 * z_roe(jj) * spd_snd_roe / rho
2861  tinv(5+jj,1) = -y_roe(jj)
2862  tinv(5+jj,5+jj) = 1.0_8
2863  enddo
2864 
2865  !check=matmul(Tmat,Tinv)
2866  !do jj = 1, 8
2867  ! do i = 1, 8
2868  ! if (i==jj.and.dabs(check(i,jj)-1.0_8)>1d-11) print*,dir,i,jj,check(i,jj)
2869  ! if (i.ne.jj.and.dabs(check(i,jj))>1d-11) print*,dir,i,jj,check(i,jj)
2870  ! enddo
2871  !enddo
2872 
2873  ! ... compute the diagonal matrix lambda
2874  lambda(:,:) = 0.0_8
2875  do jj = 1, nd
2876  lambda(jj,jj) = theta
2877  end do
2878  lambda(nd+1,nd+1) = theta + spd_snd_roe * dsqrt(xi_x**2 + xi_y**2 + xi_z**2)
2879  lambda(nd+2,nd+2) = theta - spd_snd_roe * dsqrt(xi_x**2 + xi_y**2 + xi_z**2)
2880  do jj = 1, nauxvars
2881  lambda(5+jj,5+jj) = theta
2882  enddo
2883  endif
2884  endif
2885 
2886  end subroutine sat_form_roe_matrices_multispecies
2887 
2888  subroutine pderivative_roe(RHO,en_ROE,Y_ROE,u_ref,nAuxVars,dp_drho_bar,dp_de_bar,dp_drhoY_bar)
2889 
2890 ! USE ModGlobal
2891 ! USE ModDataStruct
2892 ! USE MODMPI
2893 !#ifdef HAVE_CANTERA
2894 ! USE ModModel1
2895 ! USE Cantera
2896 !#endif
2897 
2898  Implicit none
2899  Integer :: nAuxVars
2900  Real(8) :: RHO,en_ROE,dp_drho_bar,dp_de_bar
2901  Real(8) :: Y_ROE(nauxvars),dp_drhoY_bar(nauxvars),Y_ROE_2(nauxvars+1)
2902  Real(8) :: Temp_Roe
2903 ! type(t_mixt), pointer :: state
2904  integer :: printVal
2905  Integer :: j,jj,ii
2906  Real(8) :: e_i(nauxvars),Cv_total,W_total,ja,Y_N2
2907  Real(8), pointer :: W_i(:)
2908  real(8), parameter :: R_u = 8314.4621_8 !J/K/kmole
2909  real(8) :: u_ref,Residual
2910 ! temporary placeholders for porting completion
2911  real(8) :: dimTempFactor
2912  type(t_nasapolynomials) :: polyNASA(nauxvars)
2913 
2914 ! type(t_model1), pointer:: this
2915 !#ifdef HAVE_CANTERA
2916 ! Type(t_canteraInterface), pointer :: canteraInterface
2917 !#endif
2918 ! this => state%combustion%model1
2919  y_roe_2(1:nauxvars) = y_roe(:) * rho
2920  y_roe_2(nauxvars+1) = rho
2921 !#ifdef HAVE_CANTERA
2922 ! W_i => this%canteraInterface%molW
2923 ! Temp_Roe=400.0_8
2924 ! Call poly4Temperature(rho, RHO * en_Roe, Y_ROE_2, this, Temp_ROE, printVal, errorCode=state%errorCodes)
2925 !#endif
2926  temp_roe = temp_roe * dimtempfactor ! dimensional !this%dimTempFactor
2927 
2928 
2929  If (temp_roe <= 1000.0_8) then
2930  ii = 1
2931  else
2932  ii = 2
2933  endif
2934 
2935 
2936  e_i = 0.0_8
2937  y_n2 = 1.0_8
2938  w_total = 0.0_8
2939  do jj = 1, nauxvars
2940  do j = 1, 5
2941  e_i(jj) = e_i(jj) + y_roe(jj) * polynasa(jj)%cMod(j,ii) * temp_roe**j
2942  enddo
2943  e_i(jj) = e_i(jj) + y_roe(jj) * polynasa(jj)%cMod(6,ii) !dimensionless- e_i=(e_i-e_N)
2944  y_n2 = y_n2 - y_roe(jj)
2945  w_total = w_total + y_roe(jj) / w_i(jj)
2946  enddo
2947  w_total = w_total + y_n2 / w_i(nauxvars+1)
2948  w_total = 1.0_8 / w_total !dimensional
2949 
2950  cv_total = 0.0_8
2951  do jj = 1,nauxvars !dimensional
2952  do j = 1,5
2953  cv_total = cv_total + y_roe(jj) * polynasa(jj)%coeffs(j,ii) * temp_roe**(j-1) * r_u / w_i(jj)
2954  end do
2955  cv_total = cv_total - y_roe(jj) * r_u / w_i(jj) !* Temp_Roe
2956  end do
2957 
2958  do j = 1,5
2959  cv_total = cv_total + y_n2 * polynasa(nauxvars+1)%coeffs(j,ii) * temp_roe**(j-1) * r_u / w_i(nauxvars+1)
2960  end do
2961  cv_total = cv_total - y_n2 * r_u / w_i(nauxvars+1) !* Temp_ROE !dimensional
2962 
2963  ! not that R_u / (W C_v) is dimensionless
2964  ! R_u Temp_Roe /W has the same units as u_ref**2
2965 
2966  dp_drho_bar = r_u * temp_roe / w_i(nauxvars+1) / (u_ref**2)
2967  ja = 0.0_8
2968  do jj = 1, nauxvars
2969  ja= ja + y_roe(jj) * e_i(jj)
2970  enddo
2971  dp_drho_bar = dp_drho_bar+ r_u/ (w_total*cv_total)* ja
2972 
2973  dp_de_bar = rho * r_u / (w_total * cv_total)
2974 
2975  dp_drhoy_bar = 0.0_8
2976  do jj = 1, nauxvars
2977  dp_drhoy_bar(jj) = r_u * temp_roe * (1.0_8/w_i(jj)-1.0_8/w_i(nauxvars+1))/(u_ref**2) &
2978  - r_u / (w_total * cv_total) * e_i(jj)
2979  enddo
2980 
2981  end subroutine pderivative_roe
2982 
2983  subroutine poly4temperature(rho, rhoe, Qaux, numSpecies, dimTempFactor,Tout, printVal, iterOut_out, errorOut_out, errorCode)
2984 
2985 ! USE ModGlobal
2986 ! USE ModDataStruct
2987 ! use PolynomialRoots
2988 
2989  real(8) :: Qaux(:), Tguess, hGuess, coeSum(7,2), Tout
2990  real(8) :: coe(6), Tout4, rho, rhoe, error, Tout2, Tout3
2991  real(8) :: toler, Tnext
2992  real(8) :: tLo, tHi, hLo, hHi, tTrans, hTrans, dHdT
2993  complex*16 :: sols(4)
2994  logical :: success, printMessage, useNewton
2995  integer :: printVal
2996  integer :: errorCode(4)
2997  integer, intent(out), optional :: iterOut_out
2998  integer :: iterOut
2999  real(8), intent(out), optional :: errorOut_out
3000  real(8) :: errorOut
3001 ! type(t_model1), pointer, INTENT(IN) :: this
3002  integer :: n,i, k, M
3003 ! temporary changes for porting [mtc]
3004  type(t_nasapolynomials), pointer :: poly(:)
3005  integer :: n2Index, numSpecies
3006  real(8) :: dimTempFactor
3007  integer iter
3008  real(8) :: m_transTemperature, m_minTemperature, m_maxTemperature
3009  !...Tguess & Tout are nondimensional!
3010 
3011  ! temporary for porting only [mtc]
3012  ALLOCATE(poly(numspecies+1))
3013 
3014  errorout = 0d0
3015  iterout = 0
3016  if(present(iterout_out)) iterout_out = iterout
3017  if(present(errorout_out)) errorout_out = errorout
3018 
3019  if(printval > 0) then
3020  printmessage = .true.
3021  else
3022  printmessage = .false.
3023  end if
3024 !
3025 ! poly => this%polynomials
3026 !
3027  m = ubound(qaux,1)
3028  coesum = 0d0
3029 ! Qaux(this%iN2) = rho
3030  qaux(n2index) = rho
3031 ! do i = 1, this%numspecies+1
3032  do i = 1, numspecies+1
3033 ! if(i <= M .and. i /= this%iN2) then
3034  if(i <= m .and. i /= n2index) then
3035  coesum = coesum + poly(i)%cMod*qaux(i)
3036  !coeSum = coeSum + poly(i)%cH*Qaux(i)
3037  endif
3038  enddo
3039  ! add in contribution of N2
3040  !coeSum = coeSum + poly(this%iN2)%cH*rho
3041 ! coeSum = coeSum + poly(this%iN2)%cMod*rho
3042  coesum = coesum + poly(n2index)%cMod*rho
3043 
3044  toler = 1.e-9
3045  success = .false.
3046  iter = 0
3047  error = 0.0_8
3048  errorout = error
3049  iterout = iter
3050  ttrans = m_transtemperature
3051 
3052  ! Find the enthalpy at the transition
3053  tout = m_transtemperature
3054  tout2 = tout*tout
3055  tout3 = tout2*tout
3056  tout4 = tout2*tout2
3057  coe = coesum(1:6,1)
3058  htrans = sum(coe(1:6)*[tout,tout2,tout3,tout4,tout3*tout2,1d0])
3059 
3060  !write(*,*) "rhoe, hTrans"
3061  !write(*,*) rhoe, hTrans
3062 
3063 
3064  ! Evaluate the enthalpy to see which way to go
3065  if (rhoe < htrans) then
3066  ! Evaluate the lower bound
3067  tout = m_mintemperature
3068  tout2 = tout*tout
3069  tout3 = tout2*tout
3070  tout4 = tout2*tout2
3071  coe = coesum(1:6,1)
3072  hlo = sum(coe(1:6)*[tout,tout2,tout3,tout4,tout3*tout2,1d0])
3073  if (rhoe <= hlo) then
3074  ! Right at the minimum temperature, we're done
3075  ! Tout = Tout/this%dimTempFactor
3076  tout = tout/dimtempfactor
3077  success = .true.
3078  if (rhoe < hlo) then
3079  ! Below the maximum temperature, Clamp to min and issue a warning to the user
3080  if(printmessage) then
3081  printval = printval - 1
3082  write(*,*) "Warning: Exceeded minimum temperature in poly4."
3083  call poly4errorprint(rho, rhoe, qaux, tout, 0d0, 0)
3084  endif
3085  if (errorcode(2) >= 0) then
3086  errorcode(2) = errorcode(2) + 1
3087  else
3088  errorcode(2) = errorcode(2) - 1
3089  endif
3090  endif
3091  return
3092  endif
3093  tlo = tout
3094  thi = ttrans
3095  hhi = htrans
3096  else if (rhoe > htrans) then
3097  ! Evaluate the upper bound
3098  tout = m_maxtemperature
3099  tout2 = tout*tout
3100  tout3 = tout2*tout
3101  tout4 = tout2*tout2
3102  coe = coesum(1:6,2)
3103  hhi = sum(coe(1:6)*[tout,tout2,tout3,tout4,tout3*tout2,1d0])
3104  if (rhoe >= hhi) then
3105  ! Right at the maximum temperature, we're done
3106 ! Tout = Tout/this%dimTempFactor
3107  tout = tout/dimtempfactor
3108  if (rhoe > hhi) then
3109  ! Above the maximum temperature, Clamp to max and issue a warning to the user
3110  if(printmessage) then
3111  printval = printval - 1
3112  write(*,*) "Warning: Exceeded maximum temperature in poly4."
3113  call poly4errorprint(rho, rhoe, qaux, tout, 0d0, 0)
3114  success=.true.
3115  endif
3116  if (errorcode(3) >= 0) then
3117  errorcode(3) = errorcode(3) + 1
3118  else
3119  errorcode(3) = errorcode(3) - 1
3120  endif
3121  endif
3122  return
3123  endif
3124  tlo = ttrans
3125  thi = tout
3126  hlo = htrans
3127  else
3128  ! Right at the transition temperature, we're done
3129  tout = ttrans
3130  success = .true.
3131  endif
3132 
3133  ! Safe to start the Newton iteration now, make the initial guess based on the bracketing
3134  ! temperature/enthalpies calculated above
3135  tguess = tlo + (rhoe-hlo)*(thi-tlo)/(hhi-hlo)
3136 
3137  do i = 1,30
3138  !linearization
3139  tout = tguess
3140  tout2 = tout*tout
3141  tout3 = tout2*tout
3142  tout4 = tout2*tout2
3143  hguess = sum(coe(1:6)*[tout,tout2,tout3,tout4,tout3*tout2,1d0])
3144  dhdt = sum(coe(1:5)*[1d0,2d0*tout,3d0*tout2,4d0*tout3,5d0*tout4])
3145  error = (rhoe - hguess)/dhdt
3146  tguess = tguess + error
3147 
3148  ! Don't go beyond tLo or tHi
3149  if (tguess < tlo) tguess = tlo
3150  if (tguess > thi) tguess = thi
3151 
3152  if(abs(error/tguess) < toler ) then
3153  tout = tguess
3154  iterout=i
3155  errorout = abs(error/tout)
3156  success = .true.
3157  EXIT
3158  endif
3159  enddo
3160 
3161  ! Newton iteration failed, this should not happen. Warn the user
3162  if(.not. success .or. tout < m_mintemperature .or. tout > m_maxtemperature) then
3163  if(printmessage) then
3164  printval = printval - 1
3165  write(*,*) "Warning: Failed to find temperature in poly4 with Newton iteration."
3166  call poly4errorprint(rho, rhoe, qaux, tout, errorout, iterout)
3167  endif
3168  if (errorcode(1) >= 0) then
3169  errorcode(1) = errorcode(1) + 1
3170  else
3171  errorcode(1) = errorcode(1) - 1
3172  endif
3173  endif
3174 
3175 ! Tout = Tout / this%dimTempFactor
3176  tout = tout / dimtempfactor
3177  if(present(iterout_out)) iterout_out = iterout
3178  if(present(errorout_out)) errorout_out = errorout
3179 
3180  return
3181  end subroutine poly4temperature
3182 
3183  subroutine poly4errorprint(rho, rhoe, Qaux, Tout, error, iter)
3184 
3185  Implicit None
3186  real(8) :: rho, rhoe, Tout, Qaux(:)
3187  real(8) :: error
3188  integer :: iter
3189 
3190  write(*,*) "rho, rhoe, Tout, iter, error"
3191  write(*,*) rho, rhoe, tout, iter, error
3192  write(*,*) "Qaux"
3193  write(*,*) qaux
3194 
3195  end subroutine poly4errorprint
3196 
3197 END MODULE satutil
satutil::sat_form_roe_matrices_multispecies
subroutine sat_form_roe_matrices_multispecies(u_in, u_out, p_in, p_out, MT1, JAC, u_ref, ND, dir, nAuxVars, gamma, tmat, tinv, lambda)
Definition: SATUtil.sav.f90:2538
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::sat_form_roe_matrices_test
subroutine sat_form_roe_matrices_test(ND, u_in, u_out, gamma, norm, metricData, tmat, tinv, lambda)
Definition: SATUtil.sav.f90:1990
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
satutil::sat_form_roe_matrices_modified
subroutine sat_form_roe_matrices_modified(u_in, u_out, gamma, norm, MT1, dir, ND, tmat, tinv, lambda)
Definition: SATUtil.sav.f90:2273
satutil::poly4temperature
subroutine poly4temperature(rho, rhoe, Qaux, numSpecies, dimTempFactor, Tout, printVal, iterOut_out, errorOut_out, errorCode)
Definition: SATUtil.sav.f90:2984
satutil::pderivative_roe
subroutine pderivative_roe(RHO, en_ROE, Y_ROE, u_ref, nAuxVars, dp_drho_bar, dp_de_bar, dp_drhoY_bar)
Definition: SATUtil.sav.f90:2889
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
a
void size_t int size_t int size_t int int int int double int int double double *void size_t int size_t int int int int int double int size_t size_t size_t double double *void size_t int size_t int size_t size_t int double int double double *void size_t size_t size_t double * a
Definition: OperatorKernels.H:29
grid::rectilinear
integer(kind=4), parameter rectilinear
Definition: Grid.f90:7
satutil::t_nasapolynomials
Definition: SATUtil.sav.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
error
Definition: EulerRHS.H:33
satutil::sat_form_roe_matrices
subroutine sat_form_roe_matrices(ND, u_in, u_out, gamma, norm, tmat, tinv, lambda)
Definition: SATUtil.f90:1673
satutil::poly4errorprint
subroutine poly4errorprint(rho, rhoe, Qaux, Tout, error, iter)
Definition: SATUtil.sav.f90:3184