PlasCom2  1.0
XPACC Multi-physics simluation application
PC2Euler.f90
Go to the documentation of this file.
1 MODULE euler
2 
3  USE operators
4 
5  IMPLICIT NONE
6 
7 CONTAINS
8 
9  SUBROUTINE uniformrhs( numDim, gridSizes, numPoints, opInterval, gridMetric, &
10  numStencils, numStencilValues, stencilSizes, stencilStarts, &
11  stencilOffsets, stencilWeights, stencilID, numTransport, &
12  numPointsApply,applyPoints, &
13  rhoBuffer,rhoVBuffer,rhoEBuffer,transportBuffer,velHat, &
14  pressureBuffer, rhsBuffer)
15 
16  IMPLICIT NONE
17 
18  INTEGER(KIND=4), INTENT(IN) :: numDim, numStencils, numStencilValues, numComponents, numTransport
19  INTEGER(KIND=8), INTENT(IN) :: gridSizes(numdim), numPoints, numPointsApply
20  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim), applyPoints(numpointsapply)
21  INTEGER(KIND=4), INTENT(IN) :: stencilSizes(numstencils),stencilStarts(numstencils)
22  INTEGER(KIND=4), INTENT(IN) :: stencilOffsets(numvalues)
23  INTEGER(KIND=4), INTENT(IN) :: stencilID(numpoints)
24  REAL(KIND=8), INTENT(IN) :: stencilWeights(numvalues)
25  REAL(KIND=8), INTENT(IN) :: velHat(numdim*numpoints)
26  REAL(KIND=8), INTENT(IN) :: rhoBuffer(numpoints)
27  REAL(KIND=8), INTENT(IN) :: rhoVBuffer(numdim*numpoints)
28  REAL(KIND=8), INTENT(IN) :: rhoEBuffer(numpoints)
29  REAL(KIND=8), INTENT(IN) :: transportBuffer(numtransport*numpoints)
30  REAL(KIND=8), INTENT(IN) :: pressureBuffer(numpoints)
31  REAL(KIND=8), INTENT(OUT) :: rhsBuffer(numdim*numpoints)
32 
33 
34  INTEGER :: iDim, numComponents
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
37  INTEGER(KIND=8) :: iPoint2, iPoint3, pointOffset
38 
39 
40 
41  numcomponents = 1
42 
43  ! Continuity
44  DO idim = 1,numdim
45 
46  pointoffset = (idim-1)*numpoints
47 
48  DO ipoint = 1, numpointsapply
49  pointindex = applypoints(ipoint)
50  vectorpointindex = pointindex + pointoffset
51  flux(pointindex) = rhobuffer(pointindex) * velhat(vectorpointindex)
52  END DO
53 
54  CALL applyoperator(region, ng, ifirstderiv, i, flux, dflux, .false., 0)
55  do ii = 1, nc
56 #ifdef AXISYMMETRIC
57  if(nd == 2 .and. i == 1 .and. grid%XYZ(ii,1) < grid%DRmin) dflux(ii)=0d0
58 #endif
59  rhs(ii,1) = rhs(ii,1) - dflux(ii)
60  end do
61  Call save_patch_deriv(region, ng, 1, i, dflux)
62  end do ! i
63 
64 
65 #ifndef AXISYMMETRIC
66  ! ... momentum
67  do i = 1, nd
68  ip1 = i+1
69  do j = 1, nd
70  t2mapji = t2map(j,i)
71 #ifdef ENABLE_SIMD
72  !DIR$ SIMD
73 #endif
74  do ii = 1, nc
75  flux(ii) = cv(ii,ip1) * uvwhat(ii,j) + mt1(ii,t2mapji) * dv(ii,1)
76  end do
77  call apply_operator_box(region, ng, ifirstderiv, j, flux, dflux, .false., 0)
78 #ifdef ENABLE_SIMD
79  !DIR$ SIMD
80 #endif
81  do ii = 1, nc
82  rhs(ii,ip1) = rhs(ii,ip1) - dflux(ii)
83  end do
84  Call save_patch_deriv(region, ng, i+1, j, dflux)
85  end do ! j
86  end do ! i
87 #else
88 
89  ! ... momentum
90  do i = 1, nd
91  ip1 = i+1
92  do j = 1, nd
93  t2mapji = t2map(j,i)
94 #ifdef ENABLE_SIMD
95  !DIR$ SIMD
96 #endif
97  do ii = 1, nc
98  flux(ii) = cv(ii,ip1) * uvwhat(ii,j)
99  end do
100  call apply_operator_box(region, ng, ifirstderiv, j, flux, dflux, .false., 0)
101  do ii = 1, nc
102  if(nd == 2 .and. j == 1 .and. grid%XYZ(ii,1) < grid%DRmin) dflux(ii)=0d0
103  rhs(ii,ip1) = rhs(ii,ip1) - dflux(ii)
104  end do
105  Call save_patch_deriv(region, ng, i+1, j, dflux)
106  do ii = 1, nc
107  flux(ii) = mt0(ii,t2mapji) * (dv(ii,1) - 1d0/input%GamRef)
108  end do
109  call apply_operator_box(region, ng, ifirstderiv, j, flux, dflux, .false., 0)
110  do ii = 1, nc
111  dflux(ii) = dflux(ii)*jac0(ii)*invjac(ii)
112  end do
113  do ii = 1, nc
114  rhs(ii,ip1) = rhs(ii,ip1) - dflux(ii)
115  end do
116  end do ! j
117 
118  end do ! i
119 #endif
120 
121  ! ... energy
122  do i = 1, nd
123 #ifdef ENABLE_SIMD
124  !DIR$ SIMD
125 #endif
126  do ii = 1, nc
127  flux(ii) = (cv(ii,ndp2) + dv(ii,1)) * uvwhat(ii,i) - xi_tau(ii,i) * dv(ii,1)
128  end do
129  call apply_operator_box(region, ng, ifirstderiv, i, flux, dflux, .false., 0)
130  do ii = 1, nc
131 #ifdef AXISYMMETRIC
132  if(nd == 2 .and. i == 1 .and. grid%XYZ(ii,1) < grid%DRmin) dflux(ii)=0d0
133 #endif
134  rhs(ii,ndp2) = rhs(ii,ndp2) - dflux(ii)
135  end do
136  Call save_patch_deriv(region, ng, ndp2, i, dflux)
137  end do ! i
138 
139  ! ... scalar advection
140  do k = 1, nauxvars
141  do i = 1, nd
142 #ifdef ENABLE_SIMD
143  !DIR$ SIMD
144 #endif
145  do ii = 1, nc
146  flux(ii) = auxvars(ii,k) * uvwhat(ii,i)
147  end do
148  call apply_operator_box(region, ng, ifirstderiv, i, flux, dflux, .false., 0)
149  do ii = 1, nc
150 #ifdef AXISYMMETRIC
151  if(nd == 2 .and. i == 1 .and. grid%XYZ(ii,1) < grid%DRmin) dflux(ii)=0d0
152 #endif
153  rhs_auxvars(ii,k) = rhs_auxvars(ii,k) - dflux(ii)
154  end do
155  end do ! i
156 
157  end do ! k
158 
159 #ifdef AXISYMMETRIC
160  Call ns_rhs_axisymmetric_hopital(ng, nc, nauxvars, nd, ifirstderiv, region, grid, input, flux, dflux, mt0, invjac_tau, jac0, invjac, rhs, rhs_auxvars, cv, dv, auxvars, 0)
161 #endif
162 
163  ! ... deallocate
164  deallocate( uvwhat )
165 
166  ! CAA Benchmark 2, Category 1, Problem 1
167  if (.false.) then
168 #ifdef ENABLE_SIMD
169  !DIR$ SIMD
170 #endif
171  do i = 1, nc
172  source = 0.1_dp * exp(-(log(2.0_dp))*( ((xyz(i,1)-4.0_dp)**2+xyz(i,2)**2)/(0.4_dp)**2 )) &
173  * dsin(4.0_dp * twopi * time(i)) * invjac(i)
174  rhs(i,4) = rhs(i,4) + source / (gv(i,1)-1.0_dp)
175  end do
176  end if
177 
178  ! ... subtract off (violation of) metric identity
179 
180  if (usemetricidentities == true) then
181 
182  if( nd == 2 ) then
183 
184 #ifdef ENABLE_SIMD
185  !DIR$ SIMD
186 #endif
187  do ii = 1, nc
188  rhs(ii,1) = rhs(ii,1) + cv(ii,2) * ident(ii,1) &
189  + cv(ii,3) * ident(ii,2)
190  rhs(ii,2) = rhs(ii,2) + dv(ii,1) * ident(ii,1)
191  rhs(ii,3) = rhs(ii,3) + dv(ii,1) * ident(ii,2)
192  rhs(ii,4) = rhs(ii,4) + (cv(ii,4) + dv(ii,1)) * cv(ii,2)/cv(ii,1) * ident(ii,1) &
193  + (cv(ii,4) + dv(ii,1)) * cv(ii,3)/cv(ii,1) * ident(ii,2)
194  end do
195 
196  else if( nd == 3 ) then
197 
198 #ifdef ENABLE_SIMD
199  !DIR$ SIMD
200 #endif
201  do ii = 1, nc
202  rhs(ii,1) = rhs(ii,1) + cv(ii,2) * ident(ii,1) &
203  + cv(ii,3) * ident(ii,2) &
204  + cv(ii,4) * ident(ii,3)
205  rhs(ii,2) = rhs(ii,2) + dv(ii,1) * ident(ii,1)
206  rhs(ii,3) = rhs(ii,3) + dv(ii,1) * ident(ii,2)
207  rhs(ii,4) = rhs(ii,4) + dv(ii,1) * ident(ii,3)
208  rhs(ii,5) = rhs(ii,5) + (cv(ii,5) + dv(ii,1)) * cv(ii,2)/cv(ii,1) * ident(ii,1) &
209  + (cv(ii,5) + dv(ii,1)) * cv(ii,3)/cv(ii,1) * ident(ii,2) &
210  + (cv(ii,5) + dv(ii,1)) * cv(ii,4)/cv(ii,1) * ident(ii,3)
211  end do
212 
213  end if ! ND
214 
215  end if ! useMetricIdentities
216 
217  ! ... gravity
218 
219  if (froude > 0.0_dp) then
220  norm(1) = sin(gravity_angle_phi) * cos(gravity_angle_theta) * invfroude
221  norm(2) = sin(gravity_angle_phi) * sin(gravity_angle_theta) * invfroude
222 
223  if( nd == 2 ) then
224 
225 #ifdef ENABLE_SIMD
226  !DIR$ SIMD
227 #endif
228  do i = 1, nc
229  rhs(i,2) = rhs(i,2) + cv(i,1) * norm(1) * invjac(i)
230  rhs(i,3) = rhs(i,3) + cv(i,1) * norm(2) * invjac(i)
231  end do
232 
233  else if( nd == 3 ) then
234 
235  norm(3) = cos(gravity_angle_phi) * invfroude
236 #ifdef ENABLE_SIMD
237  !DIR$ SIMD
238 #endif
239  do i = 1, nc
240  rhs(i,2) = rhs(i,2) + cv(i,1) * norm(1) * invjac(i)
241  rhs(i,3) = rhs(i,3) + cv(i,1) * norm(2) * invjac(i)
242  rhs(i,4) = rhs(i,4) + cv(i,1) * norm(3) * invjac(i)
243  end do
244 
245  end if ! ND
246 
247  end if ! Froude
248 
249  return
250 
251 
252  END SUBROUTINE eulerrhs
253 
254 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
grid
Definition: Grid.f90:1