PlasCom2  1.0
XPACC Multi-physics simluation application
Euler.f90
Go to the documentation of this file.
1 MODULE euler
2 
3  USE operators
4  USE grid
5 
6  IMPLICIT NONE
7 
8 CONTAINS
9 
10  SUBROUTINE uniformrhs( &
11  numDim, gridSizes, numPoints, &
12  fullInterval, opInterval, gridMetric, &
13  numStencils, numStencilValues, stencilSizes, stencilStarts, &
14  stencilOffsets, stencilWeights, stencilID, &
15  rhoBuffer,rhoVBuffer,rhoEBuffer,velHat, &
16  pressureBuffer, rhoRHS, rhoVRHS, rhoERHS)
17 
18  IMPLICIT NONE
19 
20  INTEGER(KIND=4), INTENT(IN) :: numDim, numStencils, numStencilValues
21  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints
22  INTEGER(KIND=8), INTENT(IN) :: fullInterval(2*numdim), opInterval(2*numdim)
23  INTEGER(KIND=4), INTENT(IN) :: stencilSizes(numstencils),stencilStarts(numstencils)
24  INTEGER(KIND=4), INTENT(IN) :: stencilOffsets(numstencilvalues)
25  INTEGER(KIND=4), INTENT(IN), TARGET :: stencilID(numdim*numpoints)
26  REAL(KIND=8), INTENT(IN) :: stencilWeights(numstencilvalues)
27  REAL(KIND=8), INTENT(IN) :: gridMetric(numdim)
28  REAL(KIND=8), INTENT(IN), TARGET :: velHat(numdim*numpoints)
29  REAL(KIND=8), INTENT(IN) :: rhoBuffer(numpoints)
30  REAL(KIND=8), INTENT(IN), TARGET :: rhoVBuffer(numdim*numpoints)
31  REAL(KIND=8), INTENT(IN) :: rhoEBuffer(numpoints)
32  REAL(KIND=8), INTENT(IN) :: pressureBuffer(numpoints)
33  REAL(KIND=8), INTENT(OUT) :: rhoRHS(numpoints)
34  REAL(KIND=8), INTENT(OUT),TARGET :: rhoVRHS(numdim*numpoints)
35  REAL(KIND=8), INTENT(OUT) :: rhoERHS(numpoints)
36 
37 
38  INTEGER :: iDim, numComponents, iVelDim
39  INTEGER(KIND=8) :: iPoint,iX,iY,iZ,iStart,iEnd,jStart,jEnd,kStart,kEnd
40  INTEGER(KIND=8) :: xIndex,zIndex,yIndex,yzIndex,xSize,ySize,zSize
41  INTEGER(KIND=8) :: iPoint2, iPoint3, pointOffset, pointIndex,vectorPointIndex
42  INTEGER(KIND=8) :: velIndex, fluxIndex, velOffset
43  REAL(KIND=8) :: gridScale
44 
45  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: fluxBuffer, dFluxBuffer, scaledPressure
46  REAL(KIND=8), DIMENSION(:), POINTER :: rhoVPtr
47  REAL(KIND=8), DIMENSION(:), POINTER :: rhoVRHSPtr
48  REAL(KIND=8), DIMENSION(:), POINTER :: bufPtr
49  INTEGER(KIND=4), DIMENSION(:), POINTER :: stencilConnectivity
50 
51  ALLOCATE(fluxbuffer(numpoints))
52  ALLOCATE(dfluxbuffer(numpoints))
53  ALLOCATE(scaledpressure(numpoints))
54 
55  numcomponents = 1
56 
57  ! Continuity
58  DO idim = 1,numdim
59 
60  pointoffset = (idim-1)*numpoints
61 
62  stencilconnectivity => stencilid((pointoffset+1):(pointoffset+numpoints))
63 
64  ! Calculates fluxBuffer = rho * vHat(iDim) [full Interval]
65  bufptr => velhat((pointoffset+1):(pointoffset+numpoints))
66  CALL zxy(numdim,numpoints,gridsizes,fullinterval,rhobuffer,bufptr,fluxbuffer)
67  ! WRITE(*,*) 'FFLUXBUFFERRHO',iDim,fluxBuffer
68  ! # need OMP barrier here to fix OMP issue?
69  ! Deriv in X(n) direction of the flux: d/dX(iDim) * [fluxBuffer]
70  CALL applyoperator(numdim, gridsizes, numcomponents, numpoints, idim, opinterval, &
71  numstencils, stencilsizes, stencilstarts, numstencilvalues, &
72  stencilweights, stenciloffsets, stencilconnectivity, fluxbuffer, dfluxbuffer)
73 
74  ! Sums dFlux into RHS [only on local interval]
75  ! rhoRHS = rhoRHS - dFlux
76  CALL yaxpy(numdim,numpoints,gridsizes,opinterval,-1.0_8,dfluxbuffer,rhorhs)
77  ! Call Save_Patch_Deriv(region, ng, 1, i, dflux)
78 
79  END DO ! iDim (Continuity)
80 
81  ! ... momentum
82  DO idim = 1, numdim
83 
84  pointoffset = (idim-1)*numpoints
85  gridscale = gridmetric(idim)
86  rhovptr => rhovbuffer((pointoffset+1):(pointoffset+numpoints))
87  rhovrhsptr => rhovrhs((pointoffset+1):(pointoffset+numpoints))
88 
89  ! Gets gridscaled pressure for diagonal (scaledPressure = gridScale*pressure)
90  CALL yax(numdim,numpoints,gridsizes,fullinterval,gridscale,pressurebuffer,scaledpressure)
91 
92  DO iveldim = 1, numdim
93 
94  veloffset = (iveldim-1)*numpoints
95  bufptr => velhat((veloffset+1):(veloffset+numpoints))
96  stencilconnectivity => stencilid((veloffset+1):(veloffset+numpoints))
97 
98  IF(idim == iveldim) THEN
99  ! Diagonal term: flux = rhoV(iDim)*vHat(iDim) + scaledPressure
100  CALL zwxpy(numdim,numpoints,gridsizes,fullinterval,rhovptr,bufptr,scaledpressure,fluxbuffer)
101  ELSE
102  ! Cross terms: flux = rhoV(iDim)*vHat(iVelDim)
103  CALL zxy(numdim,numpoints,gridsizes,fullinterval,rhovptr,bufptr,fluxbuffer)
104  ENDIF
105  ! WRITE(*,*) 'FFLUXBUFFERRHOV',iDim,iVelDim,fluxBuffer
106  ! # need OMP barrier here
107 
108  CALL applyoperator(numdim, gridsizes, numcomponents, numpoints, iveldim, opinterval, &
109  numstencils, stencilsizes, stencilstarts, numstencilvalues, &
110  stencilweights, stenciloffsets, stencilconnectivity, fluxbuffer, dfluxbuffer)
111 
112  ! Sum dFlux to RHS (using Y = a*X + Y)
113  CALL yaxpy(numdim,numpoints,gridsizes,opinterval,-1.0_8,dfluxbuffer,rhovrhsptr)
114  ! Call Save_Patch_Deriv(region, ng, i+1, j, dflux)
115  END DO ! iVelDim
116  END DO ! iDim (momentum)
117 
118 
119  ! ... energy
120  DO idim = 1, numdim
121 
122  veloffset = (idim-1)*numpoints
123  bufptr => velhat((veloffset+1):(veloffset+numpoints))
124 
125  stencilconnectivity => stencilid((veloffset+1):(veloffset+numpoints))
126 
127  ! Calculate flux = vHat(iDim)*(rhoE + pressure)
128  CALL zwmxpy(numdim,numpoints,gridsizes,fullinterval,bufptr,rhoebuffer,pressurebuffer,fluxbuffer)
129 ! WRITE(*,*) 'FFLUXBUFFERRHOE',iDim,fluxBuffer
130 
131  ! # Need OMP barrier
132  CALL applyoperator(numdim, gridsizes, numcomponents, numpoints, idim, opinterval, &
133  numstencils, stencilsizes, stencilstarts, numstencilvalues, &
134  stencilweights, stenciloffsets, stencilconnectivity, fluxbuffer, dfluxbuffer)
135 
136  CALL yaxpy(numdim,numpoints,gridsizes,opinterval,-1.0_8,dfluxbuffer,rhoerhs)
137  ! Call Save_Patch_Deriv(region, ng, i+1, j, dflux)
138  END DO ! iDim (energy)
139 
140 
141  DEALLOCATE(scaledpressure)
142  DEALLOCATE(fluxbuffer)
143  DEALLOCATE(dfluxbuffer)
144 
145  END SUBROUTINE uniformrhs
146 
147 
148 
149  SUBROUTINE uniformflux( &
150  numDim, fluxDim, gridSizes, numPoints, &
151  opInterval, gridMetric, &
152  rhoBuffer,rhoVBuffer,rhoEBuffer,velHat, &
153  pressureBuffer, scaledPressure, fluxBuffer)
154 
155  IMPLICIT NONE
156 
157  INTEGER(KIND=4), INTENT(IN) :: numDim, fluxDim
158  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints
159  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
160  REAL(KIND=8), INTENT(IN) :: gridMetric(numdim)
161  REAL(KIND=8), INTENT(IN) :: rhoBuffer(numpoints)
162  REAL(KIND=8), INTENT(IN), TARGET :: rhoVBuffer(numdim*numpoints)
163  REAL(KIND=8), INTENT(IN) :: rhoEBuffer(numpoints)
164  REAL(KIND=8), INTENT(IN), TARGET :: velHat(numdim*numpoints)
165  REAL(KIND=8), INTENT(IN) :: pressureBuffer(numpoints)
166  REAL(KIND=8), INTENT(OUT) :: scaledPressure(numpoints)
167  REAL(KIND=8), INTENT(OUT),TARGET :: fluxBuffer(numpoints*(numdim+2))
168 
169  INTEGER :: iDim, numComponents, iVelDim
170  INTEGER(KIND=8) :: iPoint,iX,iY,iZ,iStart,iEnd,jStart,jEnd,kStart,kEnd
171  INTEGER(KIND=8) :: xIndex,zIndex,yIndex,yzIndex,xSize,ySize,zSize,fluxOffset
172  INTEGER(KIND=8) :: iPoint2, iPoint3, pointOffset, pointIndex,vectorPointIndex
173  INTEGER(KIND=8) :: velIndex, fluxIndex, velOffset, dimOffset
174  REAL(KIND=8) :: gridScale
175 
176  REAL(KIND=8), DIMENSION(:), POINTER :: rhoVPtr
177  REAL(KIND=8), DIMENSION(:), POINTER :: rhoVRHSPtr
178  REAL(KIND=8), DIMENSION(:), POINTER :: bufPtr
179  REAL(KIND=8), DIMENSION(:), POINTER :: velHatPtr
180  REAL(KIND=8), DIMENSION(:), POINTER :: fluxPtr
181  INTEGER(KIND=4), DIMENSION(:), POINTER :: stencilConnectivity
182 
183  numcomponents = 1
184 
185  idim = fluxdim
186  dimoffset = (idim-1)*numpoints
187  fluxoffset = 0
188  velhatptr => velhat(dimoffset+1:dimoffset+numpoints)
189  fluxptr => fluxbuffer(1:numpoints)
190  gridscale = gridmetric(idim)
191 
192 
193  ! Continuity
194  ! Calculates fluxBuffer = rho * vHat(iDim) [full Interval]
195  CALL zxy(numdim,numpoints,gridsizes,opinterval,rhobuffer,velhatptr,fluxptr)
196  ! WRITE(*,*) 'FFLUXBUFFERRHO',iDim,fluxPtr
197  ! numDim components of Momentum
198  ! Gets gridscaled pressure for diagonal (scaledPressure = gridScale*pressure)
199  CALL yax(numdim,numpoints,gridsizes,opinterval,gridscale,pressurebuffer,scaledpressure)
200 
201  DO iveldim = 1, numdim
202 
203  fluxoffset = fluxoffset + numpoints
204  veloffset = (iveldim-1)*numpoints
205  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
206  rhovptr => rhovbuffer((veloffset+1):(veloffset+numpoints))
207 
208  IF(idim == iveldim) THEN
209  ! Diagonal term: flux = rhoV(iDim)*vHat(iDim) + scaledPressure
210  CALL zwxpy(numdim,numpoints,gridsizes,opinterval,rhovptr,velhatptr,scaledpressure,fluxptr)
211  ELSE
212  ! Cross terms: flux = rhoV(iVelDim)*vHat(iDim)
213  CALL zxy(numdim,numpoints,gridsizes,opinterval,rhovptr,velhatptr,fluxptr)
214  ENDIF
215  ! WRITE(*,*) 'FFLUXBUFFERRHOV',iVelDim,iDim,fluxPtr
216  END DO
217 
218  ! Energy
219  fluxoffset = fluxoffset + numpoints
220  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
221 
222  ! Calculate flux = vHat(iDim)*(rhoE + pressure)
223  CALL zwmxpy(numdim,numpoints,gridsizes,opinterval,velhatptr,rhoebuffer,pressurebuffer,fluxptr)
224  !WRITE(*,*) 'FFLUXBUFFERRHOE Euler',iDim,fluxPtr
225 
226  END SUBROUTINE uniformflux
227 
228  SUBROUTINE flux1d( &
229  numDim, numPoints, gridSizes, opInterval, &
230  fluxDir, gridType, gridMetric, &
231  rhoBuffer,rhoVBuffer,rhoEBuffer,velHat, &
232  pressureBuffer, fluxBuffer)
233 
234  IMPLICIT NONE
235 
236  INTEGER(KIND=4), INTENT(IN) :: numDim, fluxDir, gridType
237  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints
238  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
239  REAL(KIND=8), INTENT(IN), TARGET :: gridMetric(numdim*numdim*numpoints)
240  REAL(KIND=8), INTENT(IN) :: rhoBuffer(numpoints)
241  REAL(KIND=8), INTENT(IN), TARGET :: rhoVBuffer(numdim*numpoints)
242  REAL(KIND=8), INTENT(IN) :: rhoEBuffer(numpoints)
243  REAL(KIND=8), INTENT(IN) :: velHat(numpoints)
244  REAL(KIND=8), INTENT(IN) :: pressureBuffer(numpoints)
245  REAL(KIND=8), INTENT(OUT),TARGET :: fluxBuffer(numpoints*(numdim+2))
246 
247  INTEGER :: iDim, numComponents, iVelDim
248  INTEGER(KIND=8) :: iPoint,iX,iY,iZ,iStart,iEnd,jStart,jEnd,kStart,kEnd
249  INTEGER(KIND=8) :: xIndex,zIndex,yIndex,yzIndex,xSize,ySize,zSize,fluxOffset
250  INTEGER(KIND=8) :: iPoint2, iPoint3, pointOffset, pointIndex,vectorPointIndex
251  INTEGER(KIND=8) :: velIndex, fluxIndex, velOffset, dimOffset, metricOffset
252  REAL(KIND=8) :: gridScale
253 
254  REAL(KIND=8), DIMENSION(:), POINTER :: rhoVPtr
255  REAL(KIND=8), DIMENSION(:), POINTER :: rhoVRHSPtr
256  REAL(KIND=8), DIMENSION(:), POINTER :: fluxPtr
257  REAL(KIND=8), DIMENSION(:), POINTER :: metricPtr
258  INTEGER(KIND=4), DIMENSION(:), POINTER :: stencilConnectivity
259 
260 
261  fluxoffset = 0
262  fluxptr => fluxbuffer(1:numpoints)
263 
264  ! Continuity
265  ! Calculates fluxBuffer = rho * vHat(iDim)
266  CALL zxy(numdim,numpoints,gridsizes,opinterval,rhobuffer,velhat,fluxptr)
267 
268  ! numDim components of Momentum
269  IF(gridtype < rectilinear) THEN
270 
271  gridscale = gridmetric(fluxdir)
272 
273  DO iveldim = 1, numdim
274 
275  fluxoffset = fluxoffset + numpoints
276  veloffset = (iveldim-1)*numpoints
277  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
278  rhovptr => rhovbuffer((veloffset+1):(veloffset+numpoints))
279 
280  IF(fluxdir == iveldim) THEN
281  ! gridMetric(iDim) is single constant over all points (0 for other dimensions)
282  ! Diagonal term: flux = rhoV(iDim)*vHat(iDim) + gridMetric(iDim)*pressure
283  CALL zawpxy(numdim,numpoints,gridsizes,opinterval,gridscale,pressurebuffer,&
284  rhovptr,velhat,fluxptr)
285  ELSE
286  ! Cross terms: flux = rhoV(iVelDim)*vHat(iDim)
287  CALL zxy(numdim,numpoints,gridsizes,opinterval,rhovptr,velhat,fluxptr)
288  ENDIF
289  END DO
290 
291  ELSE IF(gridtype < curvilinear) THEN
292 
293  metricoffset = (fluxdir-1)*numpoints
294  metricptr => gridmetric(metricoffset+1:metricoffset+numpoints)
295 
296  DO iveldim = 1, numdim
297 
298  fluxoffset = fluxoffset + numpoints
299  veloffset = (iveldim-1)*numpoints
300  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
301  rhovptr => rhovbuffer((veloffset+1):(veloffset+numpoints))
302 
303  IF(fluxdir == iveldim) THEN
304  ! gridMetric(iDim) is nodal scalar (0 for other dimensions)
305  ! Diagonal term: flux = rhoV(iDim)*vHat(iDim) + gridMetric(iDim)*pressure
306  CALL zvwpxy(numdim,numpoints,gridsizes,opinterval,metricptr,pressurebuffer,rhovptr,velhat,fluxptr)
307  ELSE
308  ! Cross terms: flux = rhoV(iVelDim)*vHat(iDim)
309  CALL zxy(numdim,numpoints,gridsizes,opinterval,rhovptr,velhat,fluxptr)
310  ENDIF
311 
312  END DO
313 
314  ELSE ! gridMetric(iDim) is a nodal numDim-vector
315 
316  metricoffset = (fluxdir-1)*numdim*numpoints
317 
318  DO iveldim = 1, numdim
319 
320  fluxoffset = fluxoffset + numpoints
321  veloffset = (iveldim-1)*numpoints
322  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
323  rhovptr => rhovbuffer((veloffset+1):(veloffset+numpoints))
324  metricptr => gridmetric(metricoffset+veloffset+1:metricoffset+veloffset+numpoints)
325 
326  ! flux = rhoV(iDim)*vHat(iDim) + p*gridMetric(iDim)
327  CALL zvwpxy(numdim,numpoints,gridsizes,opinterval,metricptr,pressurebuffer,rhovptr,velhat,fluxptr)
328 
329  END DO
330 
331  ENDIF
332 
333  ! Energy
334  fluxoffset = fluxoffset + numpoints
335  fluxptr => fluxbuffer(fluxoffset+1:fluxoffset+numpoints)
336 
337  ! Calculate flux = vHat(iDim)*(rhoE + pressure)
338  CALL zwmxpy(numdim,numpoints,gridsizes,opinterval,velhat,rhoebuffer,pressurebuffer,fluxptr)
339 
340  END SUBROUTINE flux1d
341 
342  SUBROUTINE uniformscalarrhs( &
343  numDim, gridSizes, numPoints, opInterval, &
344  numStencils, numStencilValues, stencilSizes, stencilStarts, &
345  stencilOffsets, stencilWeights, stencilID, &
346  numPointsApply,applyPoints,numScalar, &
347  scalarBuffer,velHat,scalarRHS)
348 
349  IMPLICIT NONE
350 
351  INTEGER(KIND=4), INTENT(IN) :: numDim, numStencils, numStencilValues, numScalar
352  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints, numPointsApply
353  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim), applyPoints(numpointsapply)
354  INTEGER(KIND=4), INTENT(IN) :: stencilSizes(numstencils),stencilStarts(numstencils)
355  INTEGER(KIND=4), INTENT(IN) :: stencilOffsets(numstencilvalues)
356  INTEGER(KIND=4), INTENT(IN), TARGET :: stencilID(numdim*numpoints)
357  REAL(KIND=8), INTENT(IN) :: stencilWeights(numstencilvalues)
358  REAL(KIND=8), INTENT(IN) :: velHat(numdim*numpoints)
359  REAL(KIND=8), INTENT(IN) :: scalarBuffer(numscalar*numpoints)
360  REAL(KIND=8), INTENT(OUT) :: scalarRHS(numscalar*numpoints)
361 
362 
363  INTEGER :: iDim, numComponents, iVelDim, iScalar
364  INTEGER(KIND=8) :: iPoint,iX,iY,iZ,iStart,iEnd,jStart,jEnd,kStart,kEnd
365  INTEGER(KIND=8) :: xIndex,zIndex,yIndex,yzIndex,xSize,ySize,zSize
366  INTEGER(KIND=8) :: iPoint2, iPoint3, pointOffset, pointIndex,scalarIndex
367  INTEGER(KIND=8) :: velIndex, fluxIndex, velOffset, scalarOffset
368 
369 
370  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: fluxBuffer, dFluxBuffer
371  INTEGER(KIND=4), DIMENSION(:), POINTER :: stencilConnectivity
372 
373  ALLOCATE(fluxbuffer(numpoints))
374  ALLOCATE(dfluxbuffer(numpoints))
375 
376  numcomponents = 1
377 
378  ! Scalar scalar
379 
380  DO iscalar = 1, numscalar
381 
382  scalaroffset = (iscalar-1)*numpoints
383 
384  DO idim = 1,numdim
385 
386  pointoffset = (idim-1)*numpoints
387 
388  stencilconnectivity => stencilid((pointoffset+1):(pointoffset+numpoints))
389 
390  DO ipoint = 1, numpointsapply
391  fluxindex = applypoints(ipoint)
392  scalarindex = fluxindex + scalaroffset
393  fluxbuffer(fluxindex) = scalarbuffer(scalarindex) * velhat(fluxindex+pointoffset)
394  END DO
395 
396  CALL applyoperator(numdim, gridsizes, numcomponents, numpoints, idim, opinterval, &
397  numstencils, stencilsizes, stencilstarts, numstencilvalues, &
398  stencilweights, stenciloffsets, stencilconnectivity, fluxbuffer, dfluxbuffer)
399 
400  DO ipoint = 1, numpointsapply
401  fluxindex = applypoints(ipoint)
402  scalarindex = fluxindex + scalaroffset
403  scalarrhs(scalarindex) = scalarrhs(scalarindex) - dfluxbuffer(fluxindex)
404  END DO
405  ! Call Save_Patch_Deriv(region, ng, 1, i, dflux)
406  END DO ! iDim
407  END DO ! iScalar
408 
409  DEALLOCATE(fluxbuffer)
410  DEALLOCATE(dfluxbuffer)
411 
412  END SUBROUTINE uniformscalarrhs
413 
414 
415  SUBROUTINE uniformscalarflux( &
416  numDim, fluxDim, gridSizes, numPoints, opInterval, &
417  numScalar,scalarBuffer,velHat,scalarFlux)
418 
419  IMPLICIT NONE
420 
421  INTEGER(KIND=4), INTENT(IN) :: numDim, fluxDim, numScalar
422  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints
423  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
424  REAL(KIND=8), INTENT(IN), TARGET :: velHat(numdim*numpoints)
425  REAL(KIND=8), INTENT(IN), TARGET :: scalarBuffer(numscalar*numpoints)
426  REAL(KIND=8), INTENT(OUT),TARGET :: scalarFlux(numscalar*numpoints)
427 
428 
429  INTEGER :: iDim, numComponents, iVelDim, iScalar
430  INTEGER(KIND=8) :: iPoint,iX,iY,iZ,iStart,iEnd,jStart,jEnd,kStart,kEnd
431  INTEGER(KIND=8) :: xIndex,zIndex,yIndex,yzIndex,xSize,ySize,zSize
432  INTEGER(KIND=8) :: iPoint2, iPoint3, pointOffset, pointIndex,scalarIndex
433  INTEGER(KIND=8) :: velIndex, fluxIndex, velOffset, scalarOffset
434 
435  REAL(KIND=8), DIMENSION(:), POINTER :: fluxPtr,scalarPtr,velHatPtr
436 
437 
438  numcomponents = 1
439 
440  ! Scalar scalar flux
441 
442  idim = fluxdim
443  pointoffset = (idim-1)*numpoints
444  velhatptr => velhat(pointoffset+1:pointoffset+numpoints)
445 
446  DO iscalar = 1, numscalar
447 
448  scalaroffset = (iscalar-1)*numpoints
449 
450  fluxptr => scalarflux(scalaroffset+1:scalaroffset+numpoints)
451  scalarptr => scalarbuffer(scalaroffset+1:scalaroffset+numpoints)
452 
453  CALL zxy(numdim,numpoints,gridsizes,opinterval,scalarptr,velhatptr,fluxptr)
454 
455  END DO ! iScalar
456 
457  END SUBROUTINE uniformscalarflux
458 
460  SUBROUTINE scalarflux1d( &
461  numDim, numPoints, gridSizes, opInterval, &
462  numScalars,scalarBuffer,velHat, fluxBuffer)
463 
464  IMPLICIT NONE
465 
466  INTEGER(KIND=4), INTENT(IN) :: numDim, numScalars
467  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints
468  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
469  REAL(KIND=8), INTENT(IN), TARGET :: scalarBuffer(numscalars*numpoints)
470  REAL(KIND=8), INTENT(IN) :: velHat(numpoints)
471  REAL(KIND=8), INTENT(OUT),TARGET :: fluxBuffer(numscalars*numpoints)
472 
473  INTEGER :: iScalar
474  INTEGER(KIND=8) :: scalarOffset
475 
476  REAL(KIND=8), DIMENSION(:), POINTER :: fluxPtr
477  REAL(KIND=8), DIMENSION(:), POINTER :: scalarPtr
478 
479 
480  scalaroffset = 0
481  fluxptr => fluxbuffer(1:numpoints)
482 
483  ! Scalar transport
484  ! Calculates fluxBuffer = scalar * vHat
485  DO iscalar = 1, numscalars
486  scalarptr => scalarbuffer(scalaroffset+1:scalaroffset+numpoints)
487  fluxptr => fluxbuffer(scalaroffset+1:scalaroffset+numpoints)
488  CALL zxy(numdim,numpoints,gridsizes,opinterval,scalarptr,velhat,fluxptr)
489  scalaroffset = scalaroffset + numpoints
490  END DO
491 
492  END SUBROUTINE scalarflux1d
493 
494 END MODULE euler
euler::uniformrhs
subroutine uniformrhs(numDim, gridSizes, numPoints, fullInterval, opInterval, gridMetric, numStencils, numStencilValues, stencilSizes, stencilStarts, stencilOffsets, stencilWeights, stencilID, rhoBuffer, rhoVBuffer, rhoEBuffer, velHat, pressureBuffer, rhoRHS, rhoVRHS, rhoERHS)
Definition: Euler.f90:17
operators::applyoperator
subroutine applyoperator(numDim, dimSizes, numComponents, numPoints, opDir, opInterval, numStencils, stencilSizes, stencilStarts, numValues, stencilWeights, stencilOffsets, stencilID, U, dU)
applyoperator applies an operator specified as a stencil set to the provided state data ...
Definition: Operators.f90:36
euler
Definition: Euler.f90:1
euler::uniformscalarrhs
subroutine uniformscalarrhs(numDim, gridSizes, numPoints, opInterval, numStencils, numStencilValues, stencilSizes, stencilStarts, stencilOffsets, stencilWeights, stencilID, numPointsApply, applyPoints, numScalar, scalarBuffer, velHat, scalarRHS)
Definition: Euler.f90:348
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
euler::scalarflux1d
subroutine scalarflux1d(numDim, numPoints, gridSizes, opInterval, numScalars, scalarBuffer, velHat, fluxBuffer)
Flux for scalar transport.
Definition: Euler.f90:463
euler::flux1d
subroutine flux1d(numDim, numPoints, gridSizes, opInterval, fluxDir, gridType, gridMetric, rhoBuffer, rhoVBuffer, rhoEBuffer, velHat, pressureBuffer, fluxBuffer)
Definition: Euler.f90:233
operators::zawpxy
subroutine zawpxy(numDim, numPoints, bufferSize, bufferInterval, a, W, X, Y, Z)
ZAWPXY point-wise operator performing Z = aW + XY.
Definition: Operators.f90:887
euler::uniformscalarflux
subroutine uniformscalarflux(numDim, fluxDim, gridSizes, numPoints, opInterval, numScalar, scalarBuffer, velHat, scalarFlux)
Definition: Euler.f90:418
grid
Definition: Grid.f90:1
grid::curvilinear
integer(kind=4), parameter curvilinear
Definition: Grid.f90:8
grid::rectilinear
integer(kind=4), parameter rectilinear
Definition: Grid.f90:7
operators::zwmxpy
subroutine zwmxpy(numDim, numPoints, bufferSize, bufferInterval, W, X, Y, Z)
ZWMXPY point-wise operator performing Z = W(X+Y) where all are vectors.
Definition: Operators.f90:997
operators::yax
subroutine yax(numDim, numPoints, bufferSize, bufferInterval, a, X, Y)
YAX point-wise operator performing Y = aX (scalar a)
Definition: Operators.f90:628
euler::uniformflux
subroutine uniformflux(numDim, fluxDim, gridSizes, numPoints, opInterval, gridMetric, rhoBuffer, rhoVBuffer, rhoEBuffer, velHat, pressureBuffer, scaledPressure, fluxBuffer)
Definition: Euler.f90:154
operators::zxy
subroutine zxy(numDim, numPoints, bufferSize, bufferInterval, X, Y, Z)
ZXY point-wise operator performing Z = XY (all vectors)
Definition: Operators.f90:679
operators::zvwpxy
subroutine zvwpxy(numDim, numPoints, bufferSize, bufferInterval, V, W, X, Y, Z)
ZVWPXY point-wise operator performing Z = VW + XY.
Definition: Operators.f90:941