PlasCom2  1.0
XPACC Multi-physics simluation application
Viscid.f90
Go to the documentation of this file.
1 MODULE viscid
2 
3  USE operators
4  USE grid
5 
6  IMPLICIT NONE
7 
8 CONTAINS
9 
10  !
11  ! Compute the viscous fluxes in 1 dimension is a uniform cartesian
12  ! coordinate system
13  !
14  ! obsolete: see StrongFlux1D, which is generalized for curvilinear coordiantes systems
15  SUBROUTINE viscidstronguniformflux( &
16  numDim, fluxDim, gridSizes, numPoints, opInterval, velocity, &
17  gridMetric, tauOneBuffer, tauTwoBuffer, tauThreeBuffer, scaledTau, &
18  heatFluxBuffer, fluxBuffer)
19 
20  IMPLICIT NONE
21 
22  INTEGER(KIND=4), INTENT(IN) :: numDim, fluxDim
23  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints
24  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
25  REAL(KIND=8), INTENT(IN) :: gridMetric(numdim)
26  REAL(KIND=8), INTENT(IN), TARGET :: velocity(numdim*numpoints)
27  REAL(KIND=8), INTENT(IN), TARGET :: tauOneBuffer(numpoints)
28  REAL(KIND=8), INTENT(IN), TARGET :: tauTwoBuffer(numpoints)
29  REAL(KIND=8), INTENT(IN), TARGET :: tauThreeBuffer(numpoints)
30  REAL(KIND=8), INTENT(OUT) :: scaledTau(numpoints)
31  REAL(KIND=8), INTENT(IN), TARGET :: heatFluxBuffer(numpoints*numdim)
32  REAL(KIND=8), INTENT(OUT),TARGET :: fluxBuffer(numpoints*(numdim+2))
33 
34  INTEGER :: iDim, numComponents, iVelDim
35  INTEGER(KIND=8) :: iPoint,iX,iY,iZ,iStart,iEnd,jStart,jEnd,kStart,kEnd
36  INTEGER(KIND=8) :: xIndex,zIndex,yIndex,yzIndex,xSize,ySize,zSize,fluxOffset
37  INTEGER(KIND=8) :: iPoint2, iPoint3, pointOffset, pointIndex,vectorPointIndex
38  INTEGER(KIND=8) :: velIndex, fluxIndex, velOffset, heatOffset, dimOffset
39  REAL(KIND=8) :: gridScale
40  REAL(KIND=8) :: minusOne
41 
42  REAL(KIND=8), DIMENSION(:), POINTER :: bufPtr
43  REAL(KIND=8), DIMENSION(:), POINTER :: velocityPtr
44  REAL(KIND=8), DIMENSION(:), POINTER :: fluxPtr
45  REAL(KIND=8), DIMENSION(:), POINTER :: tauOnePtr
46  REAL(KIND=8), DIMENSION(:), POINTER :: tauTwoPtr
47  REAL(KIND=8), DIMENSION(:), POINTER :: tauThreePtr
48  REAL(KIND=8), DIMENSION(:), POINTER :: heatFluxPtr
49 
50  numcomponents = 1
51  minusone = 1.0_8
52 
53  idim = fluxdim
54  fluxoffset = 0
55  velocityptr => velocity(1:numpoints)
56  heatfluxptr => heatfluxbuffer((fluxdim-1)*numpoints+1:numpoints+(fluxdim-1)*numpoints)
57  gridscale = gridmetric(idim)
58  fluxptr => fluxbuffer(1:numpoints)
59  tauoneptr => tauonebuffer(1:numpoints)
60  tautwoptr => tautwobuffer(1:numpoints)
61  tauthreeptr=> tauthreebuffer(1:numpoints)
62 
63  DO iveldim = 1, numdim
64 
65  !velOffset = (iVelDim-1)*numPoints
66  ! skip density
67  fluxoffset = fluxoffset + numpoints
68  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
69  ! flux = tau(iDim)
70  ! assuming tau(1,2,3) are the j components in tau_ij
71  !gridScale = minusOne*gridMetric(iDim)*gridMetric(iVelDim)
72  if(iveldim == 1) THEN
73  ! Gets grid (metric) scaled stress tensor (scaledPressure = gridScale*pressure)
74  CALL yax(numdim,numpoints,gridsizes,opinterval,gridscale,tauoneptr,fluxptr)
75  !CALL YAX(numDim,numPoints,gridSizes,opInterval,minusOne,tauOnePtr,fluxPtr)
76  ELSEIF(iveldim == 2) THEN
77  CALL yax(numdim,numpoints,gridsizes,opinterval,gridscale,tautwoptr,fluxptr)
78  !CALL YAX(numDim,numPoints,gridSizes,opInterval,minusOne,tauTwoPtr,fluxPtr)
79  ELSEIF(iveldim == 3) THEN
80  CALL yax(numdim,numpoints,gridsizes,opinterval,gridscale,tauthreeptr,fluxptr)
81  !CALL YAX(numDim,numPoints,gridSizes,opInterval,minusOne,tauThreePtr,fluxPtr)
82  ENDIF
83  !WRITE(*,*) 'FFLUXBUFFERRHOV',iVelDim,iDim,fluxPtr
84  END DO
85 
86  ! Energy
87  fluxoffset = fluxoffset + numpoints
88  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
89  ! Calculate flux = vHat(iVelDim)*(tau(iVelDim,iDim))
90  !WRITE(*,*) 'Before ZXY scaledTau',iDim,scaledTau
91  !WRITE(*,*) 'Before ZXY velocityPtr',iDim,velocityPtr
92  !WRITE(*,*) 'Before ZXY tauOneBuffer',iDim,tauOneBuffer
93 
94  ! viscous dissipation
95  CALL zxy(numdim,numpoints,gridsizes,opinterval,velocityptr,tauonebuffer,scaledtau)
96  !WRITE(*,*) 'After ZXY scaledTau',iDim,scaledTau
97  velocityptr => velocity(numpoints+1:2*numpoints)
98  CALL zwxpy(numdim,numpoints,gridsizes,opinterval,velocityptr,tautwobuffer,scaledtau,scaledtau)
99  IF(numdim == 3) THEN
100  velocityptr => velocity(2*numpoints+1:3*numpoints)
101  CALL zwxpy(numdim,numpoints,gridsizes,opinterval,velocityptr,tauthreebuffer,scaledtau,scaledtau)
102  ENDIF
103 
104  ! apply the grid metric
105  !gridScale = minusOne*gridMetric(iDim)
106  !WRITE(*,*) 'Before FFLUXBUFFERRHOE',iDim,fluxPtr(555)
107  CALL yax(numdim,numpoints,gridsizes,opinterval,gridscale,scaledtau,fluxptr)
108  !WRITE(*,*) 'After FFLUXBUFFERRHOE Viscid',iDim,fluxPtr(555)
109 
110  ! heat conduction
111  CALL yaxpy(numdim,numpoints,gridsizes,opinterval,gridscale,heatfluxptr,fluxptr)
112  !WRITE(*,*) 'After FFLUXBUFFERRHOE Conduction',iDim,fluxPtr(555)
113 
114  END SUBROUTINE viscidstronguniformflux
115 
123  SUBROUTINE strongflux1d &
124  (numdim, fluxdir, gridsizes, numpoints, opinterval, gridtype, &
125  gridmetric, taubuffer, energybuffer, fluxbuffer)
126 
127  IMPLICIT NONE
128 
129  INTEGER(KIND=4), INTENT(IN) :: numDim, fluxDir, gridType
130  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim),numPoints
131  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
132  REAL(KIND=8), INTENT(IN), TARGET :: gridMetric(numdim*numdim*numpoints)
133  REAL(KIND=8), INTENT(IN), TARGET :: tauBuffer(numpoints*numdim*(numdim+1)/2)
134  REAL(KIND=8), INTENT(IN), TARGET :: energyBuffer(numpoints*numdim)
135  REAL(KIND=8), INTENT(OUT),TARGET :: fluxBuffer(numpoints*(numdim+2))
136 
137  INTEGER :: iDim, tensorIndex, iVel
138  INTEGER(KIND=8) :: fluxOffset, tensorOffset, dirOffset, metricOffset
139  INTEGER(KIND=8) :: velIndex, fluxIndex, velOffset, heatOffset, dimOffset
140  REAL(KIND=8) :: gridScale
141  REAL(KIND=8) :: minusOne
142 
143  REAL(KIND=8), DIMENSION(:), POINTER :: bufPtr
144  REAL(KIND=8), DIMENSION(:), POINTER :: metricPtr
145  REAL(KIND=8), DIMENSION(:), POINTER :: fluxPtr
146  REAL(KIND=8), DIMENSION(:), POINTER :: tauPtr
147  REAL(KIND=8), DIMENSION(:), POINTER :: energyPtr
148 
149  INTEGER, DIMENSION(2,2) :: index2D
150  INTEGER, DIMENSION(3,3) :: index3D
151 
152  ! Map 0-based index for (2,3)-dimensional symmetric tensor
153  index2d = reshape((/ 0, 1, 1, 2 /), shape(index2d))
154  index3d = reshape((/ 0, 1, 2, 1, 3, 4, 2, 4, 5 /), shape(index3d))
155 
156  diroffset = (fluxdir-1)*numpoints
157  fluxoffset = 0
158 
159  IF(gridtype < rectilinear) THEN
160 
161  gridscale = gridmetric(fluxdir)
162 
163  ! Momentum terms
164  DO idim = 1, numdim
165 
166  IF(numdim == 2) THEN
167  tensorindex = index2d(fluxdir,idim)
168  ELSE
169  tensorindex = index3d(fluxdir,idim)
170  ENDIF
171 
172  tensoroffset = tensorindex*numpoints
173  fluxoffset = fluxoffset + numpoints
174  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
175  tauptr => taubuffer(tensoroffset+1:tensoroffset+numpoints)
176 
177  ! Grid metric is scalar and constant, this gets flux = metric * tau_(iDim)(
178  CALL yax(numdim,numpoints,gridsizes,opinterval,gridscale,tauptr,fluxptr)
179 
180  END DO
181 
182  ! Energy
183  fluxoffset = fluxoffset + numpoints
184  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
185  energyptr => energybuffer(diroffset+1:diroffset+numpoints)
186 
187  CALL yax(numdim,numpoints,gridsizes,opinterval,gridscale,energyptr,fluxptr)
188 
189  ELSE IF (gridtype == rectilinear) THEN
190 
191  metricptr => gridmetric(diroffset+1:diroffset+numpoints)
192 
193  ! Momentum terms
194  DO idim = 1, numdim
195 
196  IF(numdim == 2) THEN
197  tensorindex = index2d(fluxdir,idim)
198  ELSE
199  tensorindex = index3d(fluxdir,idim)
200  ENDIF
201 
202  tensoroffset = tensorindex*numpoints
203  fluxoffset = fluxoffset + numpoints
204  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
205  tauptr => taubuffer(tensoroffset+1:tensoroffset+numpoints)
206 
207  ! Grid metric is scalar, this gets flux = metric * tau_(iDim)(
208  CALL zxy(numdim,numpoints,gridsizes,opinterval,metricptr,tauptr,fluxptr)
209 
210  END DO
211 
212  ! Energy
213  fluxoffset = fluxoffset + numpoints
214  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
215  energyptr => energybuffer(diroffset+1:diroffset+numpoints)
216 
217  CALL zxy(numdim,numpoints,gridsizes,opinterval,metricptr,energyptr,fluxptr)
218 
219  ELSE ! must be curvilinear
220 
221  ! need row 'fluxDir' of the metric tensor
222  metricoffset = (fluxdir-1)*numdim*numpoints
223 
224  ! For curvilinear
225  ! Momentum terms flux_i = (metric_fluxDir_j * tau_j_i) i=[1:numDim]
226  DO idim = 1, numdim
227 
228  fluxoffset = fluxoffset + numpoints
229  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
230 
231  CALL assignmentxa(numdim,numpoints,gridsizes,opinterval,0.0_8,fluxptr)
232 
233  DO ivel = 1, numdim
234 
235  veloffset = (ivel-1)*numpoints+metricoffset
236  metricptr => gridmetric(veloffset+1:veloffset+numpoints)
237 
238  IF(numdim == 2) THEN
239  tensorindex = index2d(ivel,idim)
240  ELSE
241  tensorindex = index3d(ivel,idim)
242  ENDIF
243 
244  tensoroffset = tensorindex*numpoints
245  tauptr => taubuffer(tensoroffset+1:tensoroffset+numpoints)
246 
247  ! Get fluxPtr += (metric_fluxDir_iVel * tau_iVel_iDim)
248  CALL ywxpy(numdim,numpoints,gridsizes,opinterval,metricptr,tauptr,fluxptr)
249 
250  END DO
251 
252  END DO
253 
254  ! Energy flux is just (metric_fluxDir_i * Q_i), since
255  ! both have proper vector storage the DOT operator
256  ! can be used
257  fluxoffset = fluxoffset + numpoints
258  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
259  metricptr => gridmetric(metricoffset+1:metricoffset+numdim*numpoints)
260 
261  CALL zxdoty(numdim,numpoints,gridsizes,opinterval,numdim,metricptr,energybuffer,fluxptr)
262 
263  ENDIF
264 
265  END SUBROUTINE strongflux1d
266 
274  SUBROUTINE scalarflux1d &
275  (numdim, fluxdir, gridsizes, numpoints, opinterval, gridtype, &
276  gridmetric, gradscalar, fluxbuffer)
277 
278  IMPLICIT NONE
279 
280  INTEGER(KIND=4), INTENT(IN) :: numDim, fluxDir, gridType
281  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim),numPoints
282  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
283  REAL(KIND=8), INTENT(IN), TARGET :: gridMetric(numdim*numdim*numpoints)
284  REAL(KIND=8), INTENT(IN), TARGET :: gradScalar(numpoints*numdim)
285  REAL(KIND=8), INTENT(OUT) :: fluxBuffer(numpoints)
286 
287  INTEGER(KIND=8) :: dirOffset, metricOffset
288  INTEGER(KIND=8) :: velIndex, fluxIndex, velOffset, heatOffset, dimOffset
289  REAL(KIND=8) :: gridScale
290  REAL(KIND=8) :: minusOne
291 
292  REAL(KIND=8), DIMENSION(:), POINTER :: gradPtr
293  REAL(KIND=8), DIMENSION(:), POINTER :: metricPtr
294 
295  diroffset = (fluxdir-1)*numpoints
296 
297  IF(gridtype < rectilinear) THEN
298 
299  gridscale = gridmetric(fluxdir)
300  gradptr => gradscalar(diroffset+1:diroffset+numpoints)
301 
302  ! Grid metric is scalar and constant, this gets flux = metric * grad(Y)
303  CALL yax(numdim,numpoints,gridsizes,opinterval,gridscale,gradptr,fluxbuffer)
304 
305  ELSE IF (gridtype == rectilinear) THEN
306 
307  metricptr => gridmetric(diroffset+1:diroffset+numpoints)
308  gradptr => gradscalar(diroffset+1:diroffset+numpoints)
309 
310  ! Grid metric is scalar, this gets flux = metric * tau_(iDim)(
311  CALL zxy(numdim,numpoints,gridsizes,opinterval,metricptr,gradptr,fluxbuffer)
312 
313  ELSE ! must be curvilinear
314 
315  metricoffset = diroffset*numdim
316  metricptr => gridmetric(metricoffset+1:metricoffset+numdim*numpoints)
317 
318  CALL zxdoty(numdim,numpoints,gridsizes,opinterval,numdim,metricptr,gradscalar,fluxbuffer)
319 
320  ENDIF
321 
322  END SUBROUTINE scalarflux1d
323 
324 END MODULE viscid
operators::ywxpy
subroutine ywxpy(numDim, numPoints, bufferSize, bufferInterval, W, X, Y)
YWXPY point-wise operator performing Y = WX + Y, where all are vectors.
Definition: Operators.f90:835
viscid::viscidstronguniformflux
subroutine viscidstronguniformflux(numDim, fluxDim, gridSizes, numPoints, opInterval, velocity, gridMetric, tauOneBuffer, tauTwoBuffer, tauThreeBuffer, scaledTau, heatFluxBuffer, fluxBuffer)
Definition: Viscid.f90:19
operators::yaxpy
subroutine yaxpy(numDim, numPoints, bufferSize, bufferInterval, a, X, Y)
YAXPY point-wise operator performing Y = aX + Y (scalar a)
Definition: Operators.f90:424
operators::zwxpy
subroutine zwxpy(numDim, numPoints, bufferSize, bufferInterval, W, X, Y, Z)
ZWXPY point-wise operator performing Z = WX + Y, where all are vectors.
Definition: Operators.f90:782
viscid::scalarflux1d
subroutine scalarflux1d(numDim, fluxDir, gridSizes, numPoints, opInterval, gridType, gridMetric, gradScalar, fluxBuffer)
Compute the curvilinear cartesian viscous fluxes in 1 dimension.
Definition: Viscid.f90:277
operators::assignmentxa
subroutine assignmentxa(numDim, numPoints, bufferSize, bufferInterval, a, X)
ASSIGNMENTXA point-wise operator performing X = scalar a.
Definition: Operators.f90:1222
viscid::strongflux1d
subroutine strongflux1d(numDim, fluxDir, gridSizes, numPoints, opInterval, gridType, gridMetric, tauBuffer, energyBuffer, fluxBuffer)
Compute the curvilinear cartesian viscous fluxes in 1 dimension.
Definition: Viscid.f90:126
viscid
Definition: Viscid.f90:1
operators::zxdoty
subroutine zxdoty(numDim, numPoints, bufferSize, bufferInterval, numComponents, X, Y, Z)
ZXDOTY numComponents-vector inner product Z = X * Y.
Definition: Operators.f90:1050
grid
Definition: Grid.f90:1
grid::rectilinear
integer(kind=4), parameter rectilinear
Definition: Grid.f90:7
operators::yax
subroutine yax(numDim, numPoints, bufferSize, bufferInterval, a, X, Y)
YAX point-wise operator performing Y = aX (scalar a)
Definition: Operators.f90:628
operators::zxy
subroutine zxy(numDim, numPoints, bufferSize, bufferInterval, X, Y, Z)
ZXY point-wise operator performing Z = XY (all vectors)
Definition: Operators.f90:679