45 public :: opr_cpu_dudxyz, opr_cpu_opgrad, opr_cpu_cdtp, &
47 opr_cpu_convect_scalar, opr_cpu_set_convect_rst, &
52 module subroutine opr_cpu_dudxyz(du, u, dr, ds, dt, coef)
53 type(coef_t),
intent(in),
target :: coef
54 real(kind=
rp),
intent(inout), &
55 dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, coef%msh%nelv) :: du
56 real(kind=
rp),
intent(in), &
57 dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, coef%msh%nelv) :: &
59 end subroutine opr_cpu_dudxyz
61 module subroutine opr_cpu_opgrad(ux, uy, uz, u, coef, e_start, e_end)
62 type(coef_t),
intent(in) :: coef
63 integer,
intent(in) :: e_start, e_end
64 real(kind=
rp),
intent(inout) :: ux(coef%Xh%lxyz, e_end - e_start + 1)
65 real(kind=
rp),
intent(inout) :: uy(coef%Xh%lxyz, e_end - e_start + 1)
66 real(kind=
rp),
intent(inout) :: uz(coef%Xh%lxyz, e_end - e_start + 1)
67 real(kind=
rp),
intent(in) :: u(coef%Xh%lxyz, e_end - e_start + 1)
68 end subroutine opr_cpu_opgrad
70 module subroutine opr_cpu_cdtp(dtx, x, dr, ds, dt, coef, e_start, e_end)
71 type(coef_t),
intent(in) :: coef
72 integer,
intent(in) :: e_start, e_end
73 real(kind=
rp),
intent(inout) :: dtx(coef%Xh%lxyz, e_end - e_start + 1)
74 real(kind=
rp),
intent(inout) :: x(coef%Xh%lxyz, e_end - e_start + 1)
75 real(kind=
rp),
intent(in) :: dr(coef%Xh%lxyz, e_end - e_start + 1)
76 real(kind=
rp),
intent(in) :: ds(coef%Xh%lxyz, e_end - e_start + 1)
77 real(kind=
rp),
intent(in) :: dt(coef%Xh%lxyz, e_end - e_start + 1)
78 end subroutine opr_cpu_cdtp
80 module subroutine opr_cpu_conv1(du, u, vx, vy, vz, xh, &
82 type(space_t),
intent(in) :: Xh
83 type(coef_t),
intent(in) :: coef
84 integer,
intent(in) :: e_start, e_end
85 real(kind=
rp),
intent(inout) :: du(xh%lxyz, e_end - e_start + 1)
86 real(kind=
rp),
intent(in) :: &
87 u(xh%lx, xh%ly, xh%lz, e_end - e_start + 1)
88 real(kind=
rp),
intent(in) :: &
89 vx(xh%lx, xh%ly, xh%lz, e_end - e_start + 1)
90 real(kind=
rp),
intent(in) :: &
91 vy(xh%lx, xh%ly, xh%lz, e_end - e_start + 1)
92 real(kind=
rp),
intent(in) :: &
93 vz(xh%lx, xh%ly, xh%lz, e_end - e_start + 1)
94 end subroutine opr_cpu_conv1
96 module subroutine opr_cpu_convect_scalar(du, u, cr, cs, ct, xh_gll, &
97 xh_gl, coef_gll, coef_gl, gll_to_gl)
98 type(space_t),
intent(in) :: Xh_GL
99 type(space_t),
intent(in) :: Xh_GLL
100 type(coef_t),
intent(in) :: coef_GLL
101 type(coef_t),
intent(in) :: coef_GL
102 type(interpolator_t),
intent(inout) :: GLL_to_GL
103 real(kind=
rp),
intent(inout) :: &
104 du(xh_gll%lx, xh_gll%ly, xh_gll%lz, coef_gl%msh%nelv)
105 real(kind=
rp),
intent(inout) :: &
106 u(xh_gl%lx, xh_gl%lx, xh_gl%lx, coef_gl%msh%nelv)
107 real(kind=
rp),
intent(inout) :: cr(xh_gl%lxyz, coef_gl%msh%nelv)
108 real(kind=
rp),
intent(inout) :: cs(xh_gl%lxyz, coef_gl%msh%nelv)
109 real(kind=
rp),
intent(inout) :: ct(xh_gl%lxyz, coef_gl%msh%nelv)
111 end subroutine opr_cpu_convect_scalar
113 module subroutine opr_cpu_set_convect_rst(cr, cs, ct, cx, cy, cz, &
115 type(space_t),
intent(inout) :: Xh
116 type(coef_t),
intent(inout) :: coef
117 real(kind=
rp),
dimension(Xh%lxyz, coef%msh%nelv), &
118 intent(inout) :: cr, cs, ct
119 real(kind=
rp),
dimension(Xh%lxyz, coef%msh%nelv), &
120 intent(in) :: cx, cy, cz
121 end subroutine opr_cpu_set_convect_rst
126 subroutine opr_cpu_curl(w1, w2, w3, u1, u2, u3, work1, work2, c_Xh)
127 type(
coef_t),
intent(in) :: c_xh
128 real(kind=
rp),
intent(inout), &
129 dimension(c_Xh%Xh%lx, c_Xh%Xh%ly, c_Xh%Xh%lz, c_Xh%msh%nelv) :: w1
130 real(kind=
rp),
intent(inout), &
131 dimension(c_Xh%Xh%lx, c_Xh%Xh%ly, c_Xh%Xh%lz, c_Xh%msh%nelv) :: w2
132 real(kind=
rp),
intent(inout), &
133 dimension(c_Xh%Xh%lx, c_Xh%Xh%ly, c_Xh%Xh%lz, c_Xh%msh%nelv) :: w3
134 real(kind=
rp),
intent(in), &
135 dimension(c_Xh%Xh%lx, c_Xh%Xh%ly, c_Xh%Xh%lz, c_Xh%msh%nelv) :: u1
136 real(kind=
rp),
intent(in), &
137 dimension(c_Xh%Xh%lx, c_Xh%Xh%ly, c_Xh%Xh%lz, c_Xh%msh%nelv) :: u2
138 real(kind=
rp),
intent(in), &
139 dimension(c_Xh%Xh%lx, c_Xh%Xh%ly, c_Xh%Xh%lz, c_Xh%msh%nelv) :: u3
140 real(kind=
rp),
intent(inout), &
141 dimension(c_Xh%Xh%lx, c_Xh%Xh%ly, c_Xh%Xh%lz, c_Xh%msh%nelv) :: work1
142 real(kind=
rp),
intent(inout), &
143 dimension(c_Xh%Xh%lx, c_Xh%Xh%ly, c_Xh%Xh%lz, c_Xh%msh%nelv) :: work2
150 call opr_cpu_dudxyz(work1, u3, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
151 if (gdim .eq. 3)
then
152 call opr_cpu_dudxyz(work2, u2, c_xh%drdz, c_xh%dsdz, &
154 call sub3(w1, work1, work2, n)
156 call copy(w1, work1, n)
159 if (gdim .eq. 3)
then
160 call opr_cpu_dudxyz(work1, u1, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
161 call opr_cpu_dudxyz(work2, u3, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
162 call sub3(w2, work1, work2, n)
165 call opr_cpu_dudxyz(work2, u3, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
166 call sub3(w2, work1, work2, n)
169 call opr_cpu_dudxyz(work1, u2, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
170 call opr_cpu_dudxyz(work2, u1, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
171 call sub3(w3, work1, work2, n)
174 call opcolv(w1, w2, w3, c_xh%B, gdim, n)
175 if (c_xh%cyclic)
call opr_cpu_rotate_cyc_r4(w1, w2, w3, 1, c_xh)
176 call c_xh%gs_h%op(w1, n, gs_op_add)
177 call c_xh%gs_h%op(w2, n, gs_op_add)
178 call c_xh%gs_h%op(w3, n, gs_op_add)
179 if (c_xh%cyclic)
call opr_cpu_rotate_cyc_r4(w1, w2, w3, 0, c_xh)
180 call opcolv(w1, w2, w3, c_xh%Binv, gdim, n)
182 end subroutine opr_cpu_curl
184 function opr_cpu_cfl(dt, u, v, w, Xh, coef, nelv, gdim)
result(cfl)
187 integer :: nelv, gdim
189 real(kind=
rp),
dimension(Xh%lx, Xh%ly, Xh%lz, nelv) :: u, v, w
190 real(kind=
rp) :: cflr, cfls, cflt, cflm
191 real(kind=
rp) :: ur, us, ut
193 integer :: i, j, k, e
195 if (gdim .eq. 3)
then
202 ur = ( u(i,j,k,e)*coef%drdx(i,j,k,e) &
203 + v(i,j,k,e)*coef%drdy(i,j,k,e) &
204 + w(i,j,k,e)*coef%drdz(i,j,k,e) ) &
205 * coef%jacinv(i,j,k,e)
206 us = ( u(i,j,k,e)*coef%dsdx(i,j,k,e) &
207 + v(i,j,k,e)*coef%dsdy(i,j,k,e) &
208 + w(i,j,k,e)*coef%dsdz(i,j,k,e) ) &
209 * coef%jacinv(i,j,k,e)
210 ut = ( u(i,j,k,e)*coef%dtdx(i,j,k,e) &
211 + v(i,j,k,e)*coef%dtdy(i,j,k,e) &
212 + w(i,j,k,e)*coef%dtdz(i,j,k,e) ) &
213 * coef%jacinv(i,j,k,e)
215 cflr = abs(dt*ur*xh%dr_inv(i))
216 cfls = abs(dt*us*xh%ds_inv(j))
217 cflt = abs(dt*ut*xh%dt_inv(k))
219 cflm = cflr + cfls + cflt
232 ur = ( u(i,j,1,e)*coef%drdx(i,j,1,e) &
233 + v(i,j,1,e)*coef%drdy(i,j,1,e) ) * coef%jacinv(i,j,1,e)
234 us = ( u(i,j,1,e)*coef%dsdx(i,j,1,e) &
235 + v(i,j,1,e)*coef%dsdy(i,j,1,e) ) * coef%jacinv(i,j,1,e)
237 cflr = abs(dt*ur*xh%dr_inv(i))
238 cfls = abs(dt*us*xh%ds_inv(j))
248 end function opr_cpu_cfl
250 subroutine opr_cpu_lambda2(lambda2, u, v, w, coef)
251 type(
coef_t),
intent(in) :: coef
252 real(kind=
rp),
intent(inout), &
253 dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, coef%msh%nelv) ::
lambda2
254 real(kind=
rp),
intent(in), &
255 dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, coef%msh%nelv) :: u
256 real(kind=
rp),
intent(in), &
257 dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, coef%msh%nelv) :: v
258 real(kind=
rp),
intent(in), &
259 dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, coef%msh%nelv) :: w
260 real(kind=
rp) :: grad(coef%Xh%lxyz,3,3)
262 real(kind=
xp) :: eigen(3), b, c, d, q, r, theta, l2
263 real(kind=
xp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
264 real(kind=
xp) :: a11, a22, a33, a12, a13, a23
265 real(kind=
xp) :: msk1, msk2, msk3
267 do e = 1, coef%msh%nelv
268 call opr_cpu_opgrad(grad(1,1,1), grad(1,1,2), grad(1,1,3), &
269 u(1,1,1,e), coef,e,e)
270 call opr_cpu_opgrad(grad(1,2,1), grad(1,2,2), grad(1,2,3), &
271 v(1,1,1,e), coef,e,e)
272 call opr_cpu_opgrad(grad(1,3,1), grad(1,3,2), grad(1,3,3), &
273 w(1,1,1,e), coef,e,e)
275 do i = 1, coef%Xh%lxyz
281 s12 = 0.5_xp*(grad(i,1,2) + grad(i,2,1))
282 s13 = 0.5_xp*(grad(i,1,3) + grad(i,3,1))
283 s23 = 0.5_xp*(grad(i,2,3) + grad(i,3,2))
285 o12 = 0.5_xp*(grad(i,1,2) - grad(i,2,1))
286 o13 = 0.5_xp*(grad(i,1,3) - grad(i,3,1))
287 o23 = 0.5_xp*(grad(i,2,3) - grad(i,3,2))
289 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
290 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
291 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
293 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
294 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
295 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
298 b = -(a11 + a22 + a33)
299 c = -(a12*a12 + a13*a13 + a23*a23 &
300 - a11 * a22 - a11 * a33 - a22 * a33)
301 d = -(2.0_xp * a12 * a13 * a23 - a11 * a23*a23 &
302 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
305 q = (3.0_xp * c - b*b) / 9.0_xp
306 r = (9.0_xp * c * b - 27.0_xp * d - 2.0_xp * b*b*b) / 54.0_xp
307 theta = acos( r / sqrt(-q*q*q) )
309 eigen(1) = 2.0_xp * sqrt(-q) * cos(theta / 3.0_xp) - b / 3.0_xp
310 eigen(2) = 2.0_xp * sqrt(-q) * &
311 cos((theta + 2.0_xp *
pi) / 3.0_xp) - b / 3.0_xp
312 eigen(3) = 2.0_xp * sqrt(-q) * &
313 cos((theta + 4.0_xp *
pi) / 3.0_xp) - b / 3.0_xp
314 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
315 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
316 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
317 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
318 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
319 .le. eigen(2) .and. eigen(2) .le. eigen(1))
320 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
321 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
322 .le. eigen(3) .and. eigen(3) .le. eigen(1))
324 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
329 end subroutine opr_cpu_lambda2
331 subroutine opr_cpu_rotate_cyc_r1(vx, vy, vz, idir, coef)
332 real(kind=
rp),
dimension(:),
intent(inout) :: vx, vy, vz
333 integer,
intent(in) :: idir
334 type(
coef_t),
intent(in) :: coef
335 integer :: i, j, ncyc
336 real(kind=
rp) :: vnor, vtan
338 ncyc = coef%cyc_msk(0) - 1
343 if (idir .eq. 1)
then
344 vnor = vx(j) * coef%R11(i) + vy(j) * coef%R12(i)
345 vtan = -vx(j) * coef%R12(i) + vy(j) * coef%R11(i)
346 else if (idir .eq. 0)
then
347 vnor = vx(j) * coef%R11(i) - vy(j) * coef%R12(i)
348 vtan = vx(j) * coef%R12(i) + vy(j) * coef%R11(i)
355 end subroutine opr_cpu_rotate_cyc_r1
358 subroutine opr_cpu_rotate_cyc_r4(vx, vy, vz, idir, coef)
359 real(kind=
rp),
dimension(:,:,:,:),
intent(inout) :: vx, vy, vz
360 integer,
intent(in) :: idir
361 type(
coef_t),
intent(in) :: coef
362 integer :: i, j, ncyc
363 real(kind=
rp) :: vnor, vtan
365 ncyc = coef%cyc_msk(0) - 1
370 if (idir .eq. 1)
then
371 vnor = vx(j, 1, 1, 1) * coef%R11(i) + vy(j, 1, 1, 1) * coef%R12(i)
372 vtan = -vx(j, 1, 1, 1) * coef%R12(i) + vy(j, 1, 1, 1) * coef%R11(i)
373 else if (idir .eq. 0)
then
374 vnor = vx(j, 1, 1, 1) * coef%R11(i) - vy(j, 1, 1, 1) * coef%R12(i)
375 vtan = vx(j, 1, 1, 1) * coef%R12(i) + vy(j, 1, 1, 1) * coef%R11(i)
378 vx(j, 1, 1, 1) = vnor
379 vy(j, 1, 1, 1) = vtan
383 end subroutine opr_cpu_rotate_cyc_r4
Routines to interpolate between different spaces.
A simulation component that computes lambda2 The values are stored in the field registry under the na...
real(kind=rp), parameter, public pi
subroutine, public sub3(a, b, c, n)
Vector subtraction .
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
Collection of vector field operations operating on and . Note that in general the indices and ....
subroutine, public opcolv(a1, a2, a3, c, gdim, n)
integer, parameter, public xp
integer, parameter, public dp
integer, parameter, public rp
Global precision used in computations.
subroutine, public opr_cpu_curl(w1, w2, w3, u1, u2, u3, work1, work2, c_xh)
subroutine, public opr_cpu_rotate_cyc_r1(vx, vy, vz, idir, coef)
real(kind=rp) function, public opr_cpu_cfl(dt, u, v, w, xh, coef, nelv, gdim)
subroutine, public opr_cpu_lambda2(lambda2, u, v, w, coef)
subroutine, public opr_cpu_rotate_cyc_r4(vx, vy, vz, idir, coef)
Defines a function space.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Interpolation between two space::space_t.
The function space for the SEM solution fields.