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, w2, w3, n, gs_op_add)
177 if (c_xh%cyclic)
call opr_cpu_rotate_cyc_r4(w1, w2, w3, 0, c_xh)
178 call opcolv(w1, w2, w3, c_xh%Binv, gdim, n)
180 end subroutine opr_cpu_curl
182 function opr_cpu_cfl(dt, u, v, w, Xh, coef, nelv, gdim)
result(cfl)
185 integer :: nelv, gdim
187 real(kind=
rp),
dimension(Xh%lx, Xh%ly, Xh%lz, nelv) :: u, v, w
188 real(kind=
rp) :: cflr, cfls, cflt, cflm
189 real(kind=
rp) :: ur, us, ut
191 integer :: i, j, k, e
193 if (gdim .eq. 3)
then
200 ur = ( u(i,j,k,e)*coef%drdx(i,j,k,e) &
201 + v(i,j,k,e)*coef%drdy(i,j,k,e) &
202 + w(i,j,k,e)*coef%drdz(i,j,k,e) ) &
203 * coef%jacinv(i,j,k,e)
204 us = ( u(i,j,k,e)*coef%dsdx(i,j,k,e) &
205 + v(i,j,k,e)*coef%dsdy(i,j,k,e) &
206 + w(i,j,k,e)*coef%dsdz(i,j,k,e) ) &
207 * coef%jacinv(i,j,k,e)
208 ut = ( u(i,j,k,e)*coef%dtdx(i,j,k,e) &
209 + v(i,j,k,e)*coef%dtdy(i,j,k,e) &
210 + w(i,j,k,e)*coef%dtdz(i,j,k,e) ) &
211 * coef%jacinv(i,j,k,e)
213 cflr = abs(dt*ur*xh%dr_inv(i))
214 cfls = abs(dt*us*xh%ds_inv(j))
215 cflt = abs(dt*ut*xh%dt_inv(k))
217 cflm = cflr + cfls + cflt
230 ur = ( u(i,j,1,e)*coef%drdx(i,j,1,e) &
231 + v(i,j,1,e)*coef%drdy(i,j,1,e) ) * coef%jacinv(i,j,1,e)
232 us = ( u(i,j,1,e)*coef%dsdx(i,j,1,e) &
233 + v(i,j,1,e)*coef%dsdy(i,j,1,e) ) * coef%jacinv(i,j,1,e)
235 cflr = abs(dt*ur*xh%dr_inv(i))
236 cfls = abs(dt*us*xh%ds_inv(j))
246 end function opr_cpu_cfl
248 subroutine opr_cpu_lambda2(lambda2, u, v, w, coef)
249 type(
coef_t),
intent(in) :: coef
250 real(kind=
rp),
intent(inout), &
251 dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, coef%msh%nelv) ::
lambda2
252 real(kind=
rp),
intent(in), &
253 dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, coef%msh%nelv) :: u
254 real(kind=
rp),
intent(in), &
255 dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, coef%msh%nelv) :: v
256 real(kind=
rp),
intent(in), &
257 dimension(coef%Xh%lx, coef%Xh%ly, coef%Xh%lz, coef%msh%nelv) :: w
258 real(kind=
rp) :: grad(coef%Xh%lxyz,3,3)
260 real(kind=
xp) :: eigen(3), b, c, d, q, r, theta, l2
261 real(kind=
xp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
262 real(kind=
xp) :: a11, a22, a33, a12, a13, a23
263 real(kind=
xp) :: msk1, msk2, msk3
268 do e = 1, coef%msh%nelv
269 call opr_cpu_opgrad(grad(1,1,1), grad(1,1,2), grad(1,1,3), &
270 u(1,1,1,e), coef,e,e)
271 call opr_cpu_opgrad(grad(1,2,1), grad(1,2,2), grad(1,2,3), &
272 v(1,1,1,e), coef,e,e)
273 call opr_cpu_opgrad(grad(1,3,1), grad(1,3,2), grad(1,3,3), &
274 w(1,1,1,e), coef,e,e)
276 do i = 1, coef%Xh%lxyz
282 s12 = 0.5_xp*(grad(i,1,2) + grad(i,2,1))
283 s13 = 0.5_xp*(grad(i,1,3) + grad(i,3,1))
284 s23 = 0.5_xp*(grad(i,2,3) + grad(i,3,2))
286 o12 = 0.5_xp*(grad(i,1,2) - grad(i,2,1))
287 o13 = 0.5_xp*(grad(i,1,3) - grad(i,3,1))
288 o23 = 0.5_xp*(grad(i,2,3) - grad(i,3,2))
290 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
291 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
292 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
294 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
295 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
296 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
299 b = -(a11 + a22 + a33)
300 c = -(a12*a12 + a13*a13 + a23*a23 &
301 - a11 * a22 - a11 * a33 - a22 * a33)
302 d = -(2.0_xp * a12 * a13 * a23 - a11 * a23*a23 &
303 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
306 q = (3.0_xp * c - b*b) / 9.0_xp
307 r = (9.0_xp * c * b - 27.0_xp * d - 2.0_xp * b*b*b) / 54.0_xp
308 theta = acos( r / sqrt(-q*q*q) )
310 eigen(1) = 2.0_xp * sqrt(-q) * cos(theta / 3.0_xp) - b / 3.0_xp
311 eigen(2) = 2.0_xp * sqrt(-q) * &
312 cos((theta + 2.0_xp *
pi) / 3.0_xp) - b / 3.0_xp
313 eigen(3) = 2.0_xp * sqrt(-q) * &
314 cos((theta + 4.0_xp *
pi) / 3.0_xp) - b / 3.0_xp
315 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
316 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
317 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
318 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
319 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
320 .le. eigen(2) .and. eigen(2) .le. eigen(1))
321 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
322 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
323 .le. eigen(3) .and. eigen(3) .le. eigen(1))
325 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
331 end subroutine opr_cpu_lambda2
333 subroutine opr_cpu_rotate_cyc_r1(vx, vy, vz, idir, coef)
334 real(kind=
rp),
dimension(:),
intent(inout) :: vx, vy, vz
335 integer,
intent(in) :: idir
336 type(
coef_t),
intent(in) :: coef
337 integer :: i, j, ncyc
338 real(kind=
rp) :: vnor, vtan
340 ncyc = coef%cyc_msk(0) - 1
345 if (idir .eq. 1)
then
346 vnor = vx(j) * coef%R11(i) + vy(j) * coef%R12(i)
347 vtan = -vx(j) * coef%R12(i) + vy(j) * coef%R11(i)
348 else if (idir .eq. 0)
then
349 vnor = vx(j) * coef%R11(i) - vy(j) * coef%R12(i)
350 vtan = vx(j) * coef%R12(i) + vy(j) * coef%R11(i)
357 end subroutine opr_cpu_rotate_cyc_r1
360 subroutine opr_cpu_rotate_cyc_r4(vx, vy, vz, idir, coef)
361 real(kind=
rp),
dimension(:,:,:,:),
intent(inout) :: vx, vy, vz
362 integer,
intent(in) :: idir
363 type(
coef_t),
intent(in) :: coef
364 integer :: i, j, ncyc
365 real(kind=
rp) :: vnor, vtan
367 ncyc = coef%cyc_msk(0) - 1
372 if (idir .eq. 1)
then
373 vnor = vx(j, 1, 1, 1) * coef%R11(i) + vy(j, 1, 1, 1) * coef%R12(i)
374 vtan = -vx(j, 1, 1, 1) * coef%R12(i) + vy(j, 1, 1, 1) * coef%R11(i)
375 else if (idir .eq. 0)
then
376 vnor = vx(j, 1, 1, 1) * coef%R11(i) - vy(j, 1, 1, 1) * coef%R12(i)
377 vtan = vx(j, 1, 1, 1) * coef%R12(i) + vy(j, 1, 1, 1) * coef%R11(i)
380 vx(j, 1, 1, 1) = vnor
381 vy(j, 1, 1, 1) = vtan
385 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.