PlasCom2  1.0
XPACC Multi-physics simluation application
Grid.f90
Go to the documentation of this file.
1 MODULE grid
2 
3  IMPLICIT NONE
4 
5  INTEGER(KIND=4), PARAMETER :: cartesian = 0
6  INTEGER(KIND=4), PARAMETER :: unirect = 1
7  INTEGER(KIND=4), PARAMETER :: rectilinear = 2
8  INTEGER(KIND=4), PARAMETER :: curvilinear = 3
9 
10 CONTAINS
11 
12  SUBROUTINE computegridmetrics( &
13  numDim,gridSizes,numPoints, &
14  opInterval, numStencils, &
15  numStencilValues, stencilSizes, &
16  stencilStarts, stencilOffsets, &
17  stencilWeights,stencilID, &
18  gridCoordinates, gridMetricTensor, &
19  gridJacobianDeterminants)
20 
21  USE operators
22 
23  IMPLICIT NONE
24 
25  INTEGER(KIND=4), INTENT(IN) :: numDim
26  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints
27  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
28  INTEGER(KIND=4), INTENT(IN) :: numStencils,numStencilValues
29  INTEGER(KIND=4), INTENT(IN) :: stencilSizes(numstencils)
30  INTEGER(KIND=4), INTENT(IN) :: stencilStarts(numstencils)
31  INTEGER(KIND=4), INTENT(IN) :: stencilOffsets(numstencilvalues)
32  INTEGER(KIND=4), INTENT(IN), TARGET :: stencilID(numdim*numpoints)
33  REAL(KIND=8), INTENT(IN) :: stencilWeights(numstencilvalues)
34  REAL(KIND=8), INTENT(IN), TARGET :: gridCoordinates(numdim*numpoints)
35  REAL(KIND=8), INTENT(INOUT), TARGET :: gridMetricTensor(numdim*numdim*numpoints)
36  REAL(KIND=8), INTENT(INOUT), TARGET :: gridJacobianDeterminants(2*numpoints)
37 
38  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dCoords
39  REAL(KIND=8), DIMENSION(:), POINTER :: jacPtr, jacM1Ptr
40  REAL(KIND=8), DIMENSION(:), POINTER :: gridMetricTensorPtr
41  REAL(KIND=8), DIMENSION(:), POINTER :: coordPtr
42  INTEGER(KIND=4), DIMENSION(:), POINTER :: stencilConnPtr
43 
44  INTEGER(KIND=4) :: numComponents, iDim
45 
46  END SUBROUTINE computegridmetrics
47 
48  ! NOT THREADED
49  SUBROUTINE computecurvilineargridmetrics( &
50  numDim,gridSizes,numPoints, &
51  opInterval, numStencils, &
52  numStencilValues, stencilSizes, &
53  stencilStarts, stencilOffsets, &
54  stencilWeights,stencilID, &
55  gridCoordinates, gridMetricTensor, &
56  gridJacobianDeterminants)
57 
58  USE operators
59 
60  IMPLICIT NONE
61 
62  INTEGER(KIND=4), INTENT(IN) :: numDim
63  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints
64  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
65  INTEGER(KIND=4), INTENT(IN) :: numStencils,numStencilValues
66  INTEGER(KIND=4), INTENT(IN) :: stencilSizes(numstencils)
67  INTEGER(KIND=4), INTENT(IN) :: stencilStarts(numstencils)
68  INTEGER(KIND=4), INTENT(IN) :: stencilOffsets(numstencilvalues)
69  INTEGER(KIND=4), INTENT(IN), TARGET :: stencilID(numdim*numpoints)
70  REAL(KIND=8), INTENT(IN) :: stencilWeights(numstencilvalues)
71  REAL(KIND=8), INTENT(IN), TARGET :: gridCoordinates(numdim*numpoints)
72  REAL(KIND=8), INTENT(INOUT), TARGET :: gridMetricTensor(numdim*numdim*numpoints)
73  REAL(KIND=8), INTENT(INOUT), TARGET :: gridJacobianDeterminants(2*numpoints)
74 
75  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dCoords
76  REAL(KIND=8), DIMENSION(:), POINTER :: jacPtr, jacM1Ptr
77  REAL(KIND=8), DIMENSION(:), POINTER :: m1Ptr
78  REAL(KIND=8), DIMENSION(:), POINTER :: coordPtr
79  INTEGER(KIND=4), DIMENSION(:), POINTER :: stencilConnPtr
80 
81  INTEGER(KIND=4) :: numComponents, iDim
82  INTEGER(KIND=8) :: offSet
83 
84  ALLOCATE(dcoords(numdim*numpoints))
85 
86  numcomponents = 1
87  idim = 1
88  jacm1ptr => gridjacobiandeterminants(numpoints+1:2*numpoints)
89  jacptr => gridjacobiandeterminants(1:numpoints)
90 
91 
92  IF(numdim == 1) THEN
93  CALL applyoperator(numdim, gridsizes, numcomponents, numpoints, idim, opinterval, &
94  numstencils, stencilsizes, stencilstarts, numstencilvalues, &
95  stencilweights, stenciloffsets, stencilid, gridcoordinates, jacptr)
96 
97  ! Y = a/X
98  CALL yaxm1(numdim,numpoints,gridsizes,opinterval,1.0_8,jacm1ptr,jacptr)
99 
100  ! X = a
101  CALL assignmentxa(numdim,numpoints,gridsizes,opinterval,1.0_8,gridmetrictensor)
102 
103  ELSE IF(numdim == 2) THEN
104 
105 
106  idim = 2 ! d(xi)/dx = dy/d(eta)
107  coordptr => gridcoordinates(numpoints+1:numpoints+numpoints)
108  m1ptr => gridmetrictensor(1:numpoints)
109  jacptr => gridjacobiandeterminants(1:numpoints)
110  jacm1ptr => gridjacobiandeterminants(1:numpoints)
111  stencilconnptr => stencilid(numpoints+1:numpoints+numpoints)
112 
113  ! calculate
114  CALL applyoperator(numdim,gridsizes,numcomponents,numpoints,idim,opinterval, &
115  numstencils,stencilsizes,stencilstarts,numstencilvalues, &
116  stencilweights, stenciloffsets,stencilconnptr,coordptr,jacm1ptr)
117 
118  ! 1/J = 1/dy
119  CALL assignmentyx(numdim,numpoints,gridsizes,opinterval,jacm1ptr,m1ptr)
120 
121  idim = 1 ! d(eta)/dy = dx/d(xi)
122  offset = 3*numpoints
123  coordptr => gridcoordinates(1:numpoints)
124  m1ptr => gridmetrictensor(offset+1:offset+numpoints)
125  jacptr => gridjacobiandeterminants(1:numpoints)
126  jacm1ptr => gridjacobiandeterminants(1:numpoints)
127  stencilconnptr => stencilid(1:numpoints)
128 
129  ! calculate 1/dx
130  CALL applyoperator(numdim,gridsizes,numcomponents,numpoints,idim,opinterval, &
131  numstencils,stencilsizes,stencilstarts,numstencilvalues, &
132  stencilweights, stenciloffsets,stencilconnptr,coordptr,m1ptr)
133 
134  ! 1/J = 1/dy * 1/dx
135 ! CALL YXY(numDim,numPoints,gridSizes,opInterval,jacM1Ptr,m1Ptr)
136 
137  idim = 2 ! d(xi)/dy = -dx/d(eta)
138  offset = numpoints
139  coordptr => gridcoordinates(numpoints+1:numpoints+numpoints)
140  m1ptr => gridmetrictensor(offset+1:offset+numpoints)
141  jacm1ptr => gridjacobiandeterminants(1:numpoints)
142 
143  ! calculate 1/dx
144  CALL applyoperator(numdim,gridsizes,numcomponents,numpoints,idim,opinterval, &
145  numstencils,stencilsizes,stencilstarts,numstencilvalues, &
146  stencilweights, stenciloffsets,stencilconnptr,coordptr,m1ptr)
147 
148  idim = 1 ! d(eta)/dx = -dy/d(xi)
149 
150  offset = numpoints
151  coordptr => gridcoordinates(numpoints+1:numpoints+numpoints)
152  m1ptr => gridmetrictensor(offset+1:offset+numpoints)
153  jacm1ptr => gridjacobiandeterminants(1:numpoints)
154 
155  ! calculate 1/dx
156  CALL applyoperator(numdim,gridsizes,numcomponents,numpoints,idim,opinterval, &
157  numstencils,stencilsizes,stencilstarts,numstencilvalues, &
158  stencilweights, stenciloffsets,stencilconnptr,coordptr,m1ptr)
159 
160  ELSE IF(numdim == 3) THEN
161  END IF
162 
163  DEALLOCATE(dcoords)
164 
165  END SUBROUTINE computecurvilineargridmetrics
166 
167  SUBROUTINE applyuniformgridmetric( &
168  numDim, gridSizes, numPoints, &
169  opInterval, gridMetric, vBuffer,vHat)
170 
171  IMPLICIT NONE
172 
173  INTEGER(KIND=4), INTENT(IN) :: numDim
174  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints
175  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
176  REAL(KIND=8), INTENT(IN) :: gridMetric(numdim)
177  REAL(KIND=8), INTENT(IN) :: vBuffer(numdim*numpoints)
178  REAL(KIND=8), INTENT(OUT) :: vHat(numdim*numpoints)
179 
180 
181  INTEGER :: iDim
182  INTEGER(KIND=8) :: iPoint,iX,iY,iZ,iStart,iEnd,jStart,jEnd,kStart,kEnd
183  INTEGER(KIND=8) :: xIndex,zIndex,yIndex,yzIndex,xSize,ySize,zSize
184  INTEGER(KIND=8) :: iPoint2, iPoint3, nPlane, numPoints2, numPoints3
185  REAL(KIND=8) :: xMetric,yMetric,zMetric
186 
187  xsize = gridsizes(1)
188  istart = opinterval(1)
189  iend = opinterval(2)
190  xmetric = gridmetric(1)
191 
192 
193  IF( numdim == 1 ) THEN
194 
195  DO ix = istart, iend
196  vhat(ix) = xmetric*vbuffer(ix)
197  END DO
198 
199  ELSE IF( numdim == 2 ) THEN
200 
201  ymetric = gridmetric(2)
202  ysize = gridsizes(2)
203  jstart = opinterval(3)
204  jend = opinterval(4)
205 
206  DO iy = jstart,jend
207  yindex = (iy-1)*xsize
208  DO ix = istart,iend
209  ipoint = yindex + ix
210  ipoint2 = ipoint+numpoints
211  vhat(ipoint) = xmetric * vbuffer(ipoint)
212  vhat(ipoint2) = ymetric * vbuffer(ipoint2)
213  END DO
214  END DO
215 
216  ELSE IF( numdim == 3 ) THEN
217 
218  ymetric = gridmetric(2)
219  ysize = gridsizes(2)
220  jstart = opinterval(3)
221  jend = opinterval(4)
222 
223  zmetric = gridmetric(3)
224  zsize = gridsizes(3)
225  kstart = opinterval(5)
226  kend = opinterval(6)
227 
228  nplane = xsize * ysize
229  numpoints2 = 2*numpoints
230 
231  DO iz = kstart, kend
232  zindex = (iz-1)*nplane
233  DO iy = jstart,jend
234  yzindex = zindex + (iy-1)*xsize
235  DO ix = istart, iend
236  ipoint = yzindex + ix
237  ipoint2 = ipoint + numpoints
238  ipoint3 = ipoint + numpoints2
239  vhat(ipoint) = xmetric * vbuffer(ipoint)
240  vhat(ipoint2) = ymetric * vbuffer(ipoint2)
241  vhat(ipoint3) = zmetric * vbuffer(ipoint3)
242  END DO
243  END DO
244  END DO
245  END IF
246 
247  END SUBROUTINE applyuniformgridmetric
248 
249 
250  SUBROUTINE applycartesiangridmetric( &
251  numDim, gridSizes, numPoints, &
252  opInterval, gridMetric, vBuffer, vHat)
253 
254  IMPLICIT NONE
255 
256  INTEGER(KIND=4), INTENT(IN) :: numDim
257  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints
258  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
259  REAL(KIND=8), INTENT(IN) :: gridMetric(numdim*numpoints)
260  REAL(KIND=8), INTENT(IN) :: vBuffer(numdim*numpoints)
261  REAL(KIND=8), INTENT(OUT) :: vHat(numdim*numpoints)
262 
263 
264  INTEGER :: iDim
265  INTEGER(KIND=8) :: iPoint,iX,iY,iZ,iStart,iEnd,jStart,jEnd,kStart,kEnd
266  INTEGER(KIND=8) :: xIndex,zIndex,yIndex,yzIndex,xSize,ySize,zSize
267  INTEGER(KIND=8) :: iPoint2, iPoint3, nPlane, numPoints2, numPoints3
268  REAL(KIND=8) :: xMetric,yMetric,zMetric
269 
270  xsize = gridsizes(1)
271  istart = opinterval(1)
272  iend = opinterval(2)
273 
274  IF( numdim == 1 ) THEN
275 
276  DO ix = istart, iend
277  vhat(ix) = gridmetric(ix)*vbuffer(ipoint)
278  END DO
279 
280  ELSE IF( numdim == 2 ) THEN
281 
282  ysize = gridsizes(2)
283  jstart = opinterval(3)
284  jend = opinterval(4)
285 
286  DO iy = jstart,jend
287  yindex = (iy-1)*xsize
288  DO ix = istart,iend
289  ipoint = yindex + ix
290  ipoint2 = ipoint+numpoints
291  vhat(ipoint) = gridmetric(ipoint) * vbuffer(ipoint)
292  vhat(ipoint2) = gridmetric(ipoint2) * vbuffer(ipoint2)
293  END DO
294  END DO
295 
296  ELSE IF( numdim == 3 ) THEN
297 
298  ysize = gridsizes(2)
299  jstart = opinterval(3)
300  jend = opinterval(4)
301 
302  zsize = gridsizes(3)
303  kstart = opinterval(5)
304  kend = opinterval(6)
305 
306  nplane = xsize * ysize
307  numpoints2 = 2*numpoints
308 
309  DO iz = kstart, kend
310  zindex = (iz-1)*nplane
311  DO iy = jstart,jend
312  yzindex = zindex + (iy-1)*xsize
313  DO ix = istart, iend
314  ipoint = yzindex + ix
315  ipoint2 = ipoint + numpoints
316  ipoint3 = ipoint + numpoints2
317  vhat(ipoint) = gridmetric(ipoint) * vbuffer(ipoint)
318  vhat(ipoint2) = gridmetric(ipoint2) * vbuffer(ipoint2)
319  vhat(ipoint3) = gridmetric(ipoint3) * vbuffer(ipoint3)
320  END DO
321  END DO
322  END DO
323  END IF
324 
325  END SUBROUTINE applycartesiangridmetric
326 
327 END MODULE grid
operators::assignmentyx
subroutine assignmentyx(numDim, numPoints, bufferSize, bufferInterval, X, Y)
ASSIGNMENTYX point-wise operator performing Y = X.
Definition: Operators.f90:1170
grid::applyuniformgridmetric
subroutine applyuniformgridmetric(numDim, gridSizes, numPoints, opInterval, gridMetric, vBuffer, vHat)
Definition: Grid.f90:170
grid::computecurvilineargridmetrics
subroutine computecurvilineargridmetrics(numDim, gridSizes, numPoints, opInterval, numStencils, numStencilValues, stencilSizes, stencilStarts, stencilOffsets, stencilWeights, stencilID, gridCoordinates, gridMetricTensor, gridJacobianDeterminants)
Definition: Grid.f90:57
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
ix
Defines MPI-specific parallel global and program classes.
Definition: IX/include/Application.H:15
grid::unirect
integer(kind=4), parameter unirect
Definition: Grid.f90:6
operators::assignmentxa
subroutine assignmentxa(numDim, numPoints, bufferSize, bufferInterval, a, X)
ASSIGNMENTXA point-wise operator performing X = scalar a.
Definition: Operators.f90:1222
grid::computegridmetrics
subroutine computegridmetrics(numDim, gridSizes, numPoints, opInterval, numStencils, numStencilValues, stencilSizes, stencilStarts, stencilOffsets, stencilWeights, stencilID, gridCoordinates, gridMetricTensor, gridJacobianDeterminants)
Definition: Grid.f90:20
grid::applycartesiangridmetric
subroutine applycartesiangridmetric(numDim, gridSizes, numPoints, opInterval, gridMetric, vBuffer, vHat)
Definition: Grid.f90:253
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
grid::cartesian
integer(kind=4), parameter cartesian
Definition: Grid.f90:5
operators::yaxm1
subroutine yaxm1(numDim, numPoints, bufferSize, bufferInterval, a, X, Y)
YAXM1 point-wise operator performing Y = a/X (scalar a)
Definition: Operators.f90:1396