Neko 1.99.9
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
ax_helm_compr_sx.f90
Go to the documentation of this file.
1! Copyright (c) 2026, The Neko Authors
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without
5! modification, are permitted provided that the following conditions
6! are met:
7!
8! * Redistributions of source code must retain the above copyright
9! notice, this list of conditions and the following disclaimer.
10!
11! * Redistributions in binary form must reproduce the above
12! copyright notice, this list of conditions and the following
13! disclaimer in the documentation and/or other materials provided
14! with the distribution.
15!
16! * Neither the name of the authors nor the names of its
17! contributors may be used to endorse or promote products derived
18! from this software without specific prior written permission.
19!
20! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31! POSSIBILITY OF SUCH DAMAGE.
32!
36module ax_helm_compr_sx
37 use ax_helm_sx, only : ax_helm_sx_t
38 use num_types, only : rp
39 use coefs, only : coef_t
40 use space, only : space_t
41 use mesh, only : mesh_t
42 use math, only : addcol4
43 use utils, only : neko_error
44 implicit none
45 private
46
52 type, public, extends(ax_helm_sx_t) :: ax_helm_compr_sx_t
53 contains
54 procedure, pass(this) :: compute => ax_helm_compr_sx_compute
55 end type ax_helm_compr_sx_t
56
57contains
58
59 subroutine ax_helm_compr_sx_compute(this, w, u, coef, msh, Xh)
60 class(ax_helm_compr_sx_t), intent(in) :: this
61 type(mesh_t), intent(in) :: msh
62 type(space_t), intent(in) :: Xh
63 type(coef_t), intent(in) :: coef
64 real(kind=rp), intent(inout) :: w(xh%lx, xh%ly, xh%lz, msh%nelv)
65 real(kind=rp), intent(in) :: u(xh%lx, xh%ly, xh%lz, msh%nelv)
66 integer :: m
67
68 if (.not. allocated(coef%compression_inds)) then
69 call this%ax_helm_sx_t%compute(w, u, coef, msh, xh)
70 return
71 end if
72
73 if (size(coef%compression_inds) .ne. msh%nelv) then
74 call neko_error('Compressed geometric factors do not match the mesh')
75 end if
76 m = size(coef%G11_compr, 4)
77
78 select case(xh%lx)
79 case(14)
80 call sx_ax_helm_compr_lx14(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
81 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
82 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
83 coef%G13_compr, coef%G23_compr, &
84 coef%compression_inds, msh%nelv, m)
85 case(13)
86 call sx_ax_helm_compr_lx13(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
87 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
88 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
89 coef%G13_compr, coef%G23_compr, &
90 coef%compression_inds, msh%nelv, m)
91 case(12)
92 call sx_ax_helm_compr_lx12(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
93 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
94 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
95 coef%G13_compr, coef%G23_compr, &
96 coef%compression_inds, msh%nelv, m)
97 case(11)
98 call sx_ax_helm_compr_lx11(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
99 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
100 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
101 coef%G13_compr, coef%G23_compr, &
102 coef%compression_inds, msh%nelv, m)
103 case(10)
104 call sx_ax_helm_compr_lx10(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
105 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
106 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
107 coef%G13_compr, coef%G23_compr, &
108 coef%compression_inds, msh%nelv, m)
109 case(9)
110 call sx_ax_helm_compr_lx9(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
111 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
112 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
113 coef%G13_compr, coef%G23_compr, &
114 coef%compression_inds, msh%nelv, m)
115 case(8)
116 call sx_ax_helm_compr_lx8(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
117 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
118 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
119 coef%G13_compr, coef%G23_compr, &
120 coef%compression_inds, msh%nelv, m)
121 case(7)
122 call sx_ax_helm_compr_lx7(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
123 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
124 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
125 coef%G13_compr, coef%G23_compr, &
126 coef%compression_inds, msh%nelv, m)
127 case(6)
128 call sx_ax_helm_compr_lx6(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
129 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
130 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
131 coef%G13_compr, coef%G23_compr, &
132 coef%compression_inds, msh%nelv, m)
133 case(5)
134 call sx_ax_helm_compr_lx5(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
135 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
136 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
137 coef%G13_compr, coef%G23_compr, &
138 coef%compression_inds, msh%nelv, m)
139 case(4)
140 call sx_ax_helm_compr_lx4(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
141 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
142 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
143 coef%G13_compr, coef%G23_compr, &
144 coef%compression_inds, msh%nelv, m)
145 case(3)
146 call sx_ax_helm_compr_lx3(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
147 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
148 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
149 coef%G13_compr, coef%G23_compr, &
150 coef%compression_inds, msh%nelv, m)
151 case(2)
152 call sx_ax_helm_compr_lx2(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
153 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
154 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
155 coef%G13_compr, coef%G23_compr, &
156 coef%compression_inds, msh%nelv, m)
157 case default
158 call sx_ax_helm_compr_lx(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, &
159 xh%dyt, xh%dzt, coef%h1, coef%G11_compr, &
160 coef%G22_compr, coef%G33_compr, coef%G12_compr, &
161 coef%G13_compr, coef%G23_compr, &
162 coef%compression_inds, msh%nelv, m, xh%lx)
163 end select
164
165 if (coef%ifh2) call addcol4(w, coef%h2, coef%B, u, coef%dof%size())
166
167 end subroutine ax_helm_compr_sx_compute
168
169 subroutine sx_ax_helm_compr_lx(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
170 h1, G11, G22, G33, G12, G13, G23, ind, n, m, lx)
171 integer, intent(in) :: n, m, lx
172 integer, intent(in) :: ind(n)
173 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
174 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
175 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
176 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
177 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
178 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
179 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
180 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
181 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
182 real(kind=rp), intent(in) :: dx(lx, lx)
183 real(kind=rp), intent(in) :: dy(lx, lx)
184 real(kind=rp), intent(in) :: dz(lx, lx)
185 real(kind=rp), intent(in) :: dxt(lx, lx)
186 real(kind=rp), intent(in) :: dyt(lx, lx)
187 real(kind=rp), intent(in) :: dzt(lx, lx)
188 integer :: e, c, i, j, k, jj, kk
189 real(kind=rp) :: ur(lx, lx, lx, n)
190 real(kind=rp) :: us(lx, lx, lx, n)
191 real(kind=rp) :: ut(lx, lx, lx, n)
192 real(kind=rp) :: uur(lx, lx, lx, n)
193 real(kind=rp) :: uus(lx, lx, lx, n)
194 real(kind=rp) :: uut(lx, lx, lx, n)
195 real(kind=rp) :: wr, ws, wt
196
197 do i = 1, lx
198 do jj = 1, lx * lx * n
199 wr = 0d0
200 do kk = 1, lx
201 wr = wr + dx(i,kk)*u(kk,jj,1,1)
202 end do
203 ur(i,jj,1,1) = wr
204 end do
205 end do
206
207 do k = 1, lx
208 do i = 1, lx
209 do j = 1, lx
210 do e = 1, n
211 ws = 0d0
212 !NEC$ unroll_completely
213 do kk = 1, lx
214 ws = ws + dy(j,kk)*u(i,kk,k,e)
215 end do
216 us(i,j,k,e) = ws
217 end do
218 end do
219 end do
220 end do
221
222 do j = 1, lx
223 do i = 1, lx
224 do k = 1, lx
225 do e = 1, n
226 wt = 0d0
227 !NEC$ unroll_completely
228 do kk = 1, lx
229 wt = wt + dz(k,kk)*u(i,j,kk,e)
230 end do
231 ut(i,j,k,e) = wt
232 end do
233 end do
234 end do
235 end do
236
237 do i = 1, lx * lx * lx
238 do e = 1, n
239 c = ind(e)
240 uur(i,1,1,e) = h1(i,1,1,e) * &
241 ( g11(i,1,1,c) * ur(i,1,1,e) &
242 + g12(i,1,1,c) * us(i,1,1,e) &
243 + g13(i,1,1,c) * ut(i,1,1,e))
244
245 uus(i,1,1,e) = h1(i,1,1,e) * &
246 ( g22(i,1,1,c) * us(i,1,1,e) &
247 + g12(i,1,1,c) * ur(i,1,1,e) &
248 + g23(i,1,1,c) * ut(i,1,1,e) )
249
250 uut(i,1,1,e) = h1(i,1,1,e) * &
251 ( g33(i,1,1,c) * ut(i,1,1,e) &
252 + g13(i,1,1,c) * ur(i,1,1,e) &
253 + g23(i,1,1,c) * us(i,1,1,e))
254 end do
255 end do
256
257 do i = 1, lx
258 do jj = 1, lx * lx * n
259 wr = 0d0
260 do kk = 1, lx
261 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
262 end do
263 w(i,jj,1,1) = wr
264 end do
265 end do
266
267 do k = 1, lx
268 do i = 1, lx
269 do j = 1, lx
270 do e = 1, n
271 ws = 0d0
272 !NEC$ unroll_completely
273 do kk = 1, lx
274 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
275 end do
276 w(i,j,k,e) = w(i,j,k,e) + ws
277 end do
278 end do
279 end do
280 end do
281
282 do j = 1, lx
283 do i = 1, lx
284 do k = 1, lx
285 do e = 1, n
286 wt = 0d0
287 !NEC$ unroll_completely
288 do kk = 1, lx
289 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
290 end do
291 w(i,j,k,e) = w(i,j,k,e) + wt
292 end do
293 end do
294 end do
295 end do
296
297
298 end subroutine sx_ax_helm_compr_lx
299
300 subroutine sx_ax_helm_compr_lx14(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
301 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
302 integer, parameter :: lx = 14
303 integer, intent(in) :: n, m
304 integer, intent(in) :: ind(n)
305 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
306 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
307 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
308 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
309 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
310 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
311 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
312 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
313 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
314 real(kind=rp), intent(in) :: dx(lx, lx)
315 real(kind=rp), intent(in) :: dy(lx, lx)
316 real(kind=rp), intent(in) :: dz(lx, lx)
317 real(kind=rp), intent(in) :: dxt(lx, lx)
318 real(kind=rp), intent(in) :: dyt(lx, lx)
319 real(kind=rp), intent(in) :: dzt(lx, lx)
320 integer :: e, c, i, j, k, jj, kk
321 real(kind=rp) :: ur(lx, lx, lx, n)
322 real(kind=rp) :: us(lx, lx, lx, n)
323 real(kind=rp) :: ut(lx, lx, lx, n)
324 real(kind=rp) :: uur(lx, lx, lx, n)
325 real(kind=rp) :: uus(lx, lx, lx, n)
326 real(kind=rp) :: uut(lx, lx, lx, n)
327 real(kind=rp) :: wr, ws, wt
328
329 do i = 1, lx
330 do jj = 1, lx * lx * n
331 wr = 0d0
332 do kk = 1, lx
333 wr = wr + dx(i,kk)*u(kk,jj,1,1)
334 end do
335 ur(i,jj,1,1) = wr
336 end do
337 end do
338
339 do k = 1, lx
340 do i = 1, lx
341 do j = 1, lx
342 do e = 1, n
343 ws = 0d0
344 !NEC$ unroll_completely
345 do kk = 1, lx
346 ws = ws + dy(j,kk)*u(i,kk,k,e)
347 end do
348 us(i,j,k,e) = ws
349 end do
350 end do
351 end do
352 end do
353
354 do j = 1, lx
355 do i = 1, lx
356 do k = 1, lx
357 do e = 1, n
358 wt = 0d0
359 !NEC$ unroll_completely
360 do kk = 1, lx
361 wt = wt + dz(k,kk)*u(i,j,kk,e)
362 end do
363 ut(i,j,k,e) = wt
364 end do
365 end do
366 end do
367 end do
368
369 do i = 1, lx * lx * lx
370 do e = 1, n
371 c = ind(e)
372 uur(i,1,1,e) = h1(i,1,1,e) * &
373 ( g11(i,1,1,c) * ur(i,1,1,e) &
374 + g12(i,1,1,c) * us(i,1,1,e) &
375 + g13(i,1,1,c) * ut(i,1,1,e))
376
377 uus(i,1,1,e) = h1(i,1,1,e) * &
378 ( g22(i,1,1,c) * us(i,1,1,e) &
379 + g12(i,1,1,c) * ur(i,1,1,e) &
380 + g23(i,1,1,c) * ut(i,1,1,e) )
381
382 uut(i,1,1,e) = h1(i,1,1,e) * &
383 ( g33(i,1,1,c) * ut(i,1,1,e) &
384 + g13(i,1,1,c) * ur(i,1,1,e) &
385 + g23(i,1,1,c) * us(i,1,1,e))
386 end do
387 end do
388
389 do i = 1, lx
390 do jj = 1, lx * lx * n
391 wr = 0d0
392 do kk = 1, lx
393 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
394 end do
395 w(i,jj,1,1) = wr
396 end do
397 end do
398
399 do k = 1, lx
400 do i = 1, lx
401 do j = 1, lx
402 do e = 1, n
403 ws = 0d0
404 !NEC$ unroll_completely
405 do kk = 1, lx
406 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
407 end do
408 w(i,j,k,e) = w(i,j,k,e) + ws
409 end do
410 end do
411 end do
412 end do
413
414 do j = 1, lx
415 do i = 1, lx
416 do k = 1, lx
417 do e = 1, n
418 wt = 0d0
419 !NEC$ unroll_completely
420 do kk = 1, lx
421 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
422 end do
423 w(i,j,k,e) = w(i,j,k,e) + wt
424 end do
425 end do
426 end do
427 end do
428
429
430 end subroutine sx_ax_helm_compr_lx14
431
432 subroutine sx_ax_helm_compr_lx13(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
433 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
434 integer, parameter :: lx = 13
435 integer, intent(in) :: n, m
436 integer, intent(in) :: ind(n)
437 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
438 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
439 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
440 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
441 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
442 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
443 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
444 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
445 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
446 real(kind=rp), intent(in) :: dx(lx, lx)
447 real(kind=rp), intent(in) :: dy(lx, lx)
448 real(kind=rp), intent(in) :: dz(lx, lx)
449 real(kind=rp), intent(in) :: dxt(lx, lx)
450 real(kind=rp), intent(in) :: dyt(lx, lx)
451 real(kind=rp), intent(in) :: dzt(lx, lx)
452 integer :: e, c, i, j, k, jj, kk
453 real(kind=rp) :: ur(lx, lx, lx, n)
454 real(kind=rp) :: us(lx, lx, lx, n)
455 real(kind=rp) :: ut(lx, lx, lx, n)
456 real(kind=rp) :: uur(lx, lx, lx, n)
457 real(kind=rp) :: uus(lx, lx, lx, n)
458 real(kind=rp) :: uut(lx, lx, lx, n)
459 real(kind=rp) :: wr, ws, wt
460
461 do i = 1, lx
462 do jj = 1, lx * lx * n
463 wr = 0d0
464 do kk = 1, lx
465 wr = wr + dx(i,kk)*u(kk,jj,1,1)
466 end do
467 ur(i,jj,1,1) = wr
468 end do
469 end do
470
471 do k = 1, lx
472 do i = 1, lx
473 do j = 1, lx
474 do e = 1, n
475 ws = 0d0
476 !NEC$ unroll_completely
477 do kk = 1, lx
478 ws = ws + dy(j,kk)*u(i,kk,k,e)
479 end do
480 us(i,j,k,e) = ws
481 end do
482 end do
483 end do
484 end do
485
486 do j = 1, lx
487 do i = 1, lx
488 do k = 1, lx
489 do e = 1, n
490 wt = 0d0
491 !NEC$ unroll_completely
492 do kk = 1, lx
493 wt = wt + dz(k,kk)*u(i,j,kk,e)
494 end do
495 ut(i,j,k,e) = wt
496 end do
497 end do
498 end do
499 end do
500
501 do i = 1, lx * lx * lx
502 do e = 1, n
503 c = ind(e)
504 uur(i,1,1,e) = h1(i,1,1,e) * &
505 ( g11(i,1,1,c) * ur(i,1,1,e) &
506 + g12(i,1,1,c) * us(i,1,1,e) &
507 + g13(i,1,1,c) * ut(i,1,1,e))
508
509 uus(i,1,1,e) = h1(i,1,1,e) * &
510 ( g22(i,1,1,c) * us(i,1,1,e) &
511 + g12(i,1,1,c) * ur(i,1,1,e) &
512 + g23(i,1,1,c) * ut(i,1,1,e) )
513
514 uut(i,1,1,e) = h1(i,1,1,e) * &
515 ( g33(i,1,1,c) * ut(i,1,1,e) &
516 + g13(i,1,1,c) * ur(i,1,1,e) &
517 + g23(i,1,1,c) * us(i,1,1,e))
518 end do
519 end do
520
521 do i = 1, lx
522 do jj = 1, lx * lx * n
523 wr = 0d0
524 do kk = 1, lx
525 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
526 end do
527 w(i,jj,1,1) = wr
528 end do
529 end do
530
531 do k = 1, lx
532 do i = 1, lx
533 do j = 1, lx
534 do e = 1, n
535 ws = 0d0
536 !NEC$ unroll_completely
537 do kk = 1, lx
538 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
539 end do
540 w(i,j,k,e) = w(i,j,k,e) + ws
541 end do
542 end do
543 end do
544 end do
545
546 do j = 1, lx
547 do i = 1, lx
548 do k = 1, lx
549 do e = 1, n
550 wt = 0d0
551 !NEC$ unroll_completely
552 do kk = 1, lx
553 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
554 end do
555 w(i,j,k,e) = w(i,j,k,e) + wt
556 end do
557 end do
558 end do
559 end do
560
561
562 end subroutine sx_ax_helm_compr_lx13
563
564 subroutine sx_ax_helm_compr_lx12(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
565 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
566 integer, parameter :: lx = 12
567 integer, intent(in) :: n, m
568 integer, intent(in) :: ind(n)
569 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
570 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
571 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
572 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
573 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
574 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
575 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
576 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
577 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
578 real(kind=rp), intent(in) :: dx(lx, lx)
579 real(kind=rp), intent(in) :: dy(lx, lx)
580 real(kind=rp), intent(in) :: dz(lx, lx)
581 real(kind=rp), intent(in) :: dxt(lx, lx)
582 real(kind=rp), intent(in) :: dyt(lx, lx)
583 real(kind=rp), intent(in) :: dzt(lx, lx)
584 integer :: e, c, i, j, k, jj, kk
585 real(kind=rp) :: ur(lx, lx, lx, n)
586 real(kind=rp) :: us(lx, lx, lx, n)
587 real(kind=rp) :: ut(lx, lx, lx, n)
588 real(kind=rp) :: uur(lx, lx, lx, n)
589 real(kind=rp) :: uus(lx, lx, lx, n)
590 real(kind=rp) :: uut(lx, lx, lx, n)
591 real(kind=rp) :: wr, ws, wt
592
593 do i = 1, lx
594 do jj = 1, lx * lx * n
595 wr = 0d0
596 do kk = 1, lx
597 wr = wr + dx(i,kk)*u(kk,jj,1,1)
598 end do
599 ur(i,jj,1,1) = wr
600 end do
601 end do
602
603 do k = 1, lx
604 do i = 1, lx
605 do j = 1, lx
606 do e = 1, n
607 ws = 0d0
608 !NEC$ unroll_completely
609 do kk = 1, lx
610 ws = ws + dy(j,kk)*u(i,kk,k,e)
611 end do
612 us(i,j,k,e) = ws
613 end do
614 end do
615 end do
616 end do
617
618 do j = 1, lx
619 do i = 1, lx
620 do k = 1, lx
621 do e = 1, n
622 wt = 0d0
623 !NEC$ unroll_completely
624 do kk = 1, lx
625 wt = wt + dz(k,kk)*u(i,j,kk,e)
626 end do
627 ut(i,j,k,e) = wt
628 end do
629 end do
630 end do
631 end do
632
633 do i = 1, lx * lx * lx
634 do e = 1, n
635 c = ind(e)
636 uur(i,1,1,e) = h1(i,1,1,e) * &
637 ( g11(i,1,1,c) * ur(i,1,1,e) &
638 + g12(i,1,1,c) * us(i,1,1,e) &
639 + g13(i,1,1,c) * ut(i,1,1,e))
640
641 uus(i,1,1,e) = h1(i,1,1,e) * &
642 ( g22(i,1,1,c) * us(i,1,1,e) &
643 + g12(i,1,1,c) * ur(i,1,1,e) &
644 + g23(i,1,1,c) * ut(i,1,1,e) )
645
646 uut(i,1,1,e) = h1(i,1,1,e) * &
647 ( g33(i,1,1,c) * ut(i,1,1,e) &
648 + g13(i,1,1,c) * ur(i,1,1,e) &
649 + g23(i,1,1,c) * us(i,1,1,e))
650 end do
651 end do
652
653 do i = 1, lx
654 do jj = 1, lx * lx * n
655 wr = 0d0
656 do kk = 1, lx
657 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
658 end do
659 w(i,jj,1,1) = wr
660 end do
661 end do
662
663 do k = 1, lx
664 do i = 1, lx
665 do j = 1, lx
666 do e = 1, n
667 ws = 0d0
668 !NEC$ unroll_completely
669 do kk = 1, lx
670 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
671 end do
672 w(i,j,k,e) = w(i,j,k,e) + ws
673 end do
674 end do
675 end do
676 end do
677
678 do j = 1, lx
679 do i = 1, lx
680 do k = 1, lx
681 do e = 1, n
682 wt = 0d0
683 !NEC$ unroll_completely
684 do kk = 1, lx
685 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
686 end do
687 w(i,j,k,e) = w(i,j,k,e) + wt
688 end do
689 end do
690 end do
691 end do
692
693
694 end subroutine sx_ax_helm_compr_lx12
695
696 subroutine sx_ax_helm_compr_lx11(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
697 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
698 integer, parameter :: lx = 11
699 integer, intent(in) :: n, m
700 integer, intent(in) :: ind(n)
701 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
702 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
703 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
704 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
705 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
706 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
707 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
708 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
709 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
710 real(kind=rp), intent(in) :: dx(lx, lx)
711 real(kind=rp), intent(in) :: dy(lx, lx)
712 real(kind=rp), intent(in) :: dz(lx, lx)
713 real(kind=rp), intent(in) :: dxt(lx, lx)
714 real(kind=rp), intent(in) :: dyt(lx, lx)
715 real(kind=rp), intent(in) :: dzt(lx, lx)
716 integer :: e, c, i, j, k, jj, kk
717 real(kind=rp) :: ur(lx, lx, lx, n)
718 real(kind=rp) :: us(lx, lx, lx, n)
719 real(kind=rp) :: ut(lx, lx, lx, n)
720 real(kind=rp) :: uur(lx, lx, lx, n)
721 real(kind=rp) :: uus(lx, lx, lx, n)
722 real(kind=rp) :: uut(lx, lx, lx, n)
723 real(kind=rp) :: wr, ws, wt
724
725 do i = 1, lx
726 do jj = 1, lx * lx * n
727 wr = 0d0
728 do kk = 1, lx
729 wr = wr + dx(i,kk)*u(kk,jj,1,1)
730 end do
731 ur(i,jj,1,1) = wr
732 end do
733 end do
734
735 do k = 1, lx
736 do i = 1, lx
737 do j = 1, lx
738 do e = 1, n
739 ws = 0d0
740 !NEC$ unroll_completely
741 do kk = 1, lx
742 ws = ws + dy(j,kk)*u(i,kk,k,e)
743 end do
744 us(i,j,k,e) = ws
745 end do
746 end do
747 end do
748 end do
749
750 do j = 1, lx
751 do i = 1, lx
752 do k = 1, lx
753 do e = 1, n
754 wt = 0d0
755 !NEC$ unroll_completely
756 do kk = 1, lx
757 wt = wt + dz(k,kk)*u(i,j,kk,e)
758 end do
759 ut(i,j,k,e) = wt
760 end do
761 end do
762 end do
763 end do
764
765 do i = 1, lx * lx * lx
766 do e = 1, n
767 c = ind(e)
768 uur(i,1,1,e) = h1(i,1,1,e) * &
769 ( g11(i,1,1,c) * ur(i,1,1,e) &
770 + g12(i,1,1,c) * us(i,1,1,e) &
771 + g13(i,1,1,c) * ut(i,1,1,e))
772
773 uus(i,1,1,e) = h1(i,1,1,e) * &
774 ( g22(i,1,1,c) * us(i,1,1,e) &
775 + g12(i,1,1,c) * ur(i,1,1,e) &
776 + g23(i,1,1,c) * ut(i,1,1,e) )
777
778 uut(i,1,1,e) = h1(i,1,1,e) * &
779 ( g33(i,1,1,c) * ut(i,1,1,e) &
780 + g13(i,1,1,c) * ur(i,1,1,e) &
781 + g23(i,1,1,c) * us(i,1,1,e))
782 end do
783 end do
784
785 do i = 1, lx
786 do jj = 1, lx * lx * n
787 wr = 0d0
788 do kk = 1, lx
789 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
790 end do
791 w(i,jj,1,1) = wr
792 end do
793 end do
794
795 do k = 1, lx
796 do i = 1, lx
797 do j = 1, lx
798 do e = 1, n
799 ws = 0d0
800 !NEC$ unroll_completely
801 do kk = 1, lx
802 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
803 end do
804 w(i,j,k,e) = w(i,j,k,e) + ws
805 end do
806 end do
807 end do
808 end do
809
810 do j = 1, lx
811 do i = 1, lx
812 do k = 1, lx
813 do e = 1, n
814 wt = 0d0
815 !NEC$ unroll_completely
816 do kk = 1, lx
817 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
818 end do
819 w(i,j,k,e) = w(i,j,k,e) + wt
820 end do
821 end do
822 end do
823 end do
824
825
826 end subroutine sx_ax_helm_compr_lx11
827
828 subroutine sx_ax_helm_compr_lx10(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
829 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
830 integer, parameter :: lx = 10
831 integer, intent(in) :: n, m
832 integer, intent(in) :: ind(n)
833 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
834 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
835 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
836 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
837 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
838 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
839 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
840 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
841 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
842 real(kind=rp), intent(in) :: dx(lx, lx)
843 real(kind=rp), intent(in) :: dy(lx, lx)
844 real(kind=rp), intent(in) :: dz(lx, lx)
845 real(kind=rp), intent(in) :: dxt(lx, lx)
846 real(kind=rp), intent(in) :: dyt(lx, lx)
847 real(kind=rp), intent(in) :: dzt(lx, lx)
848 integer :: e, c, i, j, k, jj, kk
849 real(kind=rp) :: ur(lx, lx, lx, n)
850 real(kind=rp) :: us(lx, lx, lx, n)
851 real(kind=rp) :: ut(lx, lx, lx, n)
852 real(kind=rp) :: uur(lx, lx, lx, n)
853 real(kind=rp) :: uus(lx, lx, lx, n)
854 real(kind=rp) :: uut(lx, lx, lx, n)
855 real(kind=rp) :: wr, ws, wt
856
857 do i = 1, lx
858 do jj = 1, lx * lx * n
859 wr = 0d0
860 do kk = 1, lx
861 wr = wr + dx(i,kk)*u(kk,jj,1,1)
862 end do
863 ur(i,jj,1,1) = wr
864 end do
865 end do
866
867 do k = 1, lx
868 do i = 1, lx
869 do j = 1, lx
870 do e = 1, n
871 ws = 0d0
872 !NEC$ unroll_completely
873 do kk = 1, lx
874 ws = ws + dy(j,kk)*u(i,kk,k,e)
875 end do
876 us(i,j,k,e) = ws
877 end do
878 end do
879 end do
880 end do
881
882 do j = 1, lx
883 do i = 1, lx
884 do k = 1, lx
885 do e = 1, n
886 wt = 0d0
887 !NEC$ unroll_completely
888 do kk = 1, lx
889 wt = wt + dz(k,kk)*u(i,j,kk,e)
890 end do
891 ut(i,j,k,e) = wt
892 end do
893 end do
894 end do
895 end do
896
897 do i = 1, lx * lx * lx
898 do e = 1, n
899 c = ind(e)
900 uur(i,1,1,e) = h1(i,1,1,e) * &
901 ( g11(i,1,1,c) * ur(i,1,1,e) &
902 + g12(i,1,1,c) * us(i,1,1,e) &
903 + g13(i,1,1,c) * ut(i,1,1,e))
904
905 uus(i,1,1,e) = h1(i,1,1,e) * &
906 ( g22(i,1,1,c) * us(i,1,1,e) &
907 + g12(i,1,1,c) * ur(i,1,1,e) &
908 + g23(i,1,1,c) * ut(i,1,1,e) )
909
910 uut(i,1,1,e) = h1(i,1,1,e) * &
911 ( g33(i,1,1,c) * ut(i,1,1,e) &
912 + g13(i,1,1,c) * ur(i,1,1,e) &
913 + g23(i,1,1,c) * us(i,1,1,e))
914 end do
915 end do
916
917 do i = 1, lx
918 do jj = 1, lx * lx * n
919 wr = 0d0
920 do kk = 1, lx
921 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
922 end do
923 w(i,jj,1,1) = wr
924 end do
925 end do
926
927 do k = 1, lx
928 do i = 1, lx
929 do j = 1, lx
930 do e = 1, n
931 ws = 0d0
932 !NEC$ unroll_completely
933 do kk = 1, lx
934 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
935 end do
936 w(i,j,k,e) = w(i,j,k,e) + ws
937 end do
938 end do
939 end do
940 end do
941
942 do j = 1, lx
943 do i = 1, lx
944 do k = 1, lx
945 do e = 1, n
946 wt = 0d0
947 !NEC$ unroll_completely
948 do kk = 1, lx
949 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
950 end do
951 w(i,j,k,e) = w(i,j,k,e) + wt
952 end do
953 end do
954 end do
955 end do
956
957
958 end subroutine sx_ax_helm_compr_lx10
959
960 subroutine sx_ax_helm_compr_lx9(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
961 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
962 integer, parameter :: lx = 9
963 integer, intent(in) :: n, m
964 integer, intent(in) :: ind(n)
965 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
966 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
967 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
968 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
969 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
970 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
971 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
972 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
973 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
974 real(kind=rp), intent(in) :: dx(lx, lx)
975 real(kind=rp), intent(in) :: dy(lx, lx)
976 real(kind=rp), intent(in) :: dz(lx, lx)
977 real(kind=rp), intent(in) :: dxt(lx, lx)
978 real(kind=rp), intent(in) :: dyt(lx, lx)
979 real(kind=rp), intent(in) :: dzt(lx, lx)
980 integer :: e, c, i, j, k, jj, kk
981 real(kind=rp) :: ur(lx, lx, lx, n)
982 real(kind=rp) :: us(lx, lx, lx, n)
983 real(kind=rp) :: ut(lx, lx, lx, n)
984 real(kind=rp) :: uur(lx, lx, lx, n)
985 real(kind=rp) :: uus(lx, lx, lx, n)
986 real(kind=rp) :: uut(lx, lx, lx, n)
987 real(kind=rp) :: wr, ws, wt
988
989 do i = 1, lx
990 do jj = 1, lx * lx * n
991 wr = 0d0
992 do kk = 1, lx
993 wr = wr + dx(i,kk)*u(kk,jj,1,1)
994 end do
995 ur(i,jj,1,1) = wr
996 end do
997 end do
998
999 do k = 1, lx
1000 do i = 1, lx
1001 do j = 1, lx
1002 do e = 1, n
1003 ws = 0d0
1004 !NEC$ unroll_completely
1005 do kk = 1, lx
1006 ws = ws + dy(j,kk)*u(i,kk,k,e)
1007 end do
1008 us(i,j,k,e) = ws
1009 end do
1010 end do
1011 end do
1012 end do
1013
1014 do j = 1, lx
1015 do i = 1, lx
1016 do k = 1, lx
1017 do e = 1, n
1018 wt = 0d0
1019 !NEC$ unroll_completely
1020 do kk = 1, lx
1021 wt = wt + dz(k,kk)*u(i,j,kk,e)
1022 end do
1023 ut(i,j,k,e) = wt
1024 end do
1025 end do
1026 end do
1027 end do
1028
1029 do i = 1, lx * lx * lx
1030 do e = 1, n
1031 c = ind(e)
1032 uur(i,1,1,e) = h1(i,1,1,e) * &
1033 ( g11(i,1,1,c) * ur(i,1,1,e) &
1034 + g12(i,1,1,c) * us(i,1,1,e) &
1035 + g13(i,1,1,c) * ut(i,1,1,e))
1036
1037 uus(i,1,1,e) = h1(i,1,1,e) * &
1038 ( g22(i,1,1,c) * us(i,1,1,e) &
1039 + g12(i,1,1,c) * ur(i,1,1,e) &
1040 + g23(i,1,1,c) * ut(i,1,1,e) )
1041
1042 uut(i,1,1,e) = h1(i,1,1,e) * &
1043 ( g33(i,1,1,c) * ut(i,1,1,e) &
1044 + g13(i,1,1,c) * ur(i,1,1,e) &
1045 + g23(i,1,1,c) * us(i,1,1,e))
1046 end do
1047 end do
1048
1049 do i = 1, lx
1050 do jj = 1, lx * lx * n
1051 wr = 0d0
1052 do kk = 1, lx
1053 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
1054 end do
1055 w(i,jj,1,1) = wr
1056 end do
1057 end do
1058
1059 do k = 1, lx
1060 do i = 1, lx
1061 do j = 1, lx
1062 do e = 1, n
1063 ws = 0d0
1064 !NEC$ unroll_completely
1065 do kk = 1, lx
1066 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
1067 end do
1068 w(i,j,k,e) = w(i,j,k,e) + ws
1069 end do
1070 end do
1071 end do
1072 end do
1073
1074 do j = 1, lx
1075 do i = 1, lx
1076 do k = 1, lx
1077 do e = 1, n
1078 wt = 0d0
1079 !NEC$ unroll_completely
1080 do kk = 1, lx
1081 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
1082 end do
1083 w(i,j,k,e) = w(i,j,k,e) + wt
1084 end do
1085 end do
1086 end do
1087 end do
1088
1089
1090 end subroutine sx_ax_helm_compr_lx9
1091
1092 subroutine sx_ax_helm_compr_lx8(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1093 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
1094 integer, parameter :: lx = 8
1095 integer, intent(in) :: n, m
1096 integer, intent(in) :: ind(n)
1097 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1098 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1099 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1100 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
1101 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
1102 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
1103 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
1104 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
1105 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
1106 real(kind=rp), intent(in) :: dx(lx, lx)
1107 real(kind=rp), intent(in) :: dy(lx, lx)
1108 real(kind=rp), intent(in) :: dz(lx, lx)
1109 real(kind=rp), intent(in) :: dxt(lx, lx)
1110 real(kind=rp), intent(in) :: dyt(lx, lx)
1111 real(kind=rp), intent(in) :: dzt(lx, lx)
1112 integer :: e, c, i, j, k, jj, kk
1113 real(kind=rp) :: ur(lx, lx, lx, n)
1114 real(kind=rp) :: us(lx, lx, lx, n)
1115 real(kind=rp) :: ut(lx, lx, lx, n)
1116 real(kind=rp) :: uur(lx, lx, lx, n)
1117 real(kind=rp) :: uus(lx, lx, lx, n)
1118 real(kind=rp) :: uut(lx, lx, lx, n)
1119 real(kind=rp) :: wr, ws, wt
1120
1121 do i = 1, lx
1122 do jj = 1, lx * lx * n
1123 wr = 0d0
1124 do kk = 1, lx
1125 wr = wr + dx(i,kk)*u(kk,jj,1,1)
1126 end do
1127 ur(i,jj,1,1) = wr
1128 end do
1129 end do
1130
1131 do k = 1, lx
1132 do i = 1, lx
1133 do j = 1, lx
1134 do e = 1, n
1135 ws = 0d0
1136 !NEC$ unroll_completely
1137 do kk = 1, lx
1138 ws = ws + dy(j,kk)*u(i,kk,k,e)
1139 end do
1140 us(i,j,k,e) = ws
1141 end do
1142 end do
1143 end do
1144 end do
1145
1146 do j = 1, lx
1147 do i = 1, lx
1148 do k = 1, lx
1149 do e = 1, n
1150 wt = 0d0
1151 !NEC$ unroll_completely
1152 do kk = 1, lx
1153 wt = wt + dz(k,kk)*u(i,j,kk,e)
1154 end do
1155 ut(i,j,k,e) = wt
1156 end do
1157 end do
1158 end do
1159 end do
1160
1161 do i = 1, lx * lx * lx
1162 do e = 1, n
1163 c = ind(e)
1164 uur(i,1,1,e) = h1(i,1,1,e) * &
1165 ( g11(i,1,1,c) * ur(i,1,1,e) &
1166 + g12(i,1,1,c) * us(i,1,1,e) &
1167 + g13(i,1,1,c) * ut(i,1,1,e))
1168
1169 uus(i,1,1,e) = h1(i,1,1,e) * &
1170 ( g22(i,1,1,c) * us(i,1,1,e) &
1171 + g12(i,1,1,c) * ur(i,1,1,e) &
1172 + g23(i,1,1,c) * ut(i,1,1,e) )
1173
1174 uut(i,1,1,e) = h1(i,1,1,e) * &
1175 ( g33(i,1,1,c) * ut(i,1,1,e) &
1176 + g13(i,1,1,c) * ur(i,1,1,e) &
1177 + g23(i,1,1,c) * us(i,1,1,e))
1178 end do
1179 end do
1180
1181 do i = 1, lx
1182 do jj = 1, lx * lx * n
1183 wr = 0d0
1184 do kk = 1, lx
1185 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
1186 end do
1187 w(i,jj,1,1) = wr
1188 end do
1189 end do
1190
1191 do k = 1, lx
1192 do i = 1, lx
1193 do j = 1, lx
1194 do e = 1, n
1195 ws = 0d0
1196 !NEC$ unroll_completely
1197 do kk = 1, lx
1198 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
1199 end do
1200 w(i,j,k,e) = w(i,j,k,e) + ws
1201 end do
1202 end do
1203 end do
1204 end do
1205
1206 do j = 1, lx
1207 do i = 1, lx
1208 do k = 1, lx
1209 do e = 1, n
1210 wt = 0d0
1211 !NEC$ unroll_completely
1212 do kk = 1, lx
1213 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
1214 end do
1215 w(i,j,k,e) = w(i,j,k,e) + wt
1216 end do
1217 end do
1218 end do
1219 end do
1220
1221
1222 end subroutine sx_ax_helm_compr_lx8
1223
1224 subroutine sx_ax_helm_compr_lx7(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1225 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
1226 integer, parameter :: lx = 7
1227 integer, intent(in) :: n, m
1228 integer, intent(in) :: ind(n)
1229 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1230 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1231 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1232 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
1233 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
1234 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
1235 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
1236 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
1237 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
1238 real(kind=rp), intent(in) :: dx(lx, lx)
1239 real(kind=rp), intent(in) :: dy(lx, lx)
1240 real(kind=rp), intent(in) :: dz(lx, lx)
1241 real(kind=rp), intent(in) :: dxt(lx, lx)
1242 real(kind=rp), intent(in) :: dyt(lx, lx)
1243 real(kind=rp), intent(in) :: dzt(lx, lx)
1244 integer :: e, c, i, j, k, jj, kk
1245 real(kind=rp) :: ur(lx, lx, lx, n)
1246 real(kind=rp) :: us(lx, lx, lx, n)
1247 real(kind=rp) :: ut(lx, lx, lx, n)
1248 real(kind=rp) :: uur(lx, lx, lx, n)
1249 real(kind=rp) :: uus(lx, lx, lx, n)
1250 real(kind=rp) :: uut(lx, lx, lx, n)
1251 real(kind=rp) :: wr, ws, wt
1252
1253 do i = 1, lx
1254 do jj = 1, lx * lx * n
1255 wr = 0d0
1256 do kk = 1, lx
1257 wr = wr + dx(i,kk)*u(kk,jj,1,1)
1258 end do
1259 ur(i,jj,1,1) = wr
1260 end do
1261 end do
1262
1263 do k = 1, lx
1264 do i = 1, lx
1265 do j = 1, lx
1266 do e = 1, n
1267 ws = 0d0
1268 !NEC$ unroll_completely
1269 do kk = 1, lx
1270 ws = ws + dy(j,kk)*u(i,kk,k,e)
1271 end do
1272 us(i,j,k,e) = ws
1273 end do
1274 end do
1275 end do
1276 end do
1277
1278 do j = 1, lx
1279 do i = 1, lx
1280 do k = 1, lx
1281 do e = 1, n
1282 wt = 0d0
1283 !NEC$ unroll_completely
1284 do kk = 1, lx
1285 wt = wt + dz(k,kk)*u(i,j,kk,e)
1286 end do
1287 ut(i,j,k,e) = wt
1288 end do
1289 end do
1290 end do
1291 end do
1292
1293 do i = 1, lx * lx * lx
1294 do e = 1, n
1295 c = ind(e)
1296 uur(i,1,1,e) = h1(i,1,1,e) * &
1297 ( g11(i,1,1,c) * ur(i,1,1,e) &
1298 + g12(i,1,1,c) * us(i,1,1,e) &
1299 + g13(i,1,1,c) * ut(i,1,1,e))
1300
1301 uus(i,1,1,e) = h1(i,1,1,e) * &
1302 ( g22(i,1,1,c) * us(i,1,1,e) &
1303 + g12(i,1,1,c) * ur(i,1,1,e) &
1304 + g23(i,1,1,c) * ut(i,1,1,e) )
1305
1306 uut(i,1,1,e) = h1(i,1,1,e) * &
1307 ( g33(i,1,1,c) * ut(i,1,1,e) &
1308 + g13(i,1,1,c) * ur(i,1,1,e) &
1309 + g23(i,1,1,c) * us(i,1,1,e))
1310 end do
1311 end do
1312
1313 do i = 1, lx
1314 do jj = 1, lx * lx * n
1315 wr = 0d0
1316 do kk = 1, lx
1317 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
1318 end do
1319 w(i,jj,1,1) = wr
1320 end do
1321 end do
1322
1323 do k = 1, lx
1324 do i = 1, lx
1325 do j = 1, lx
1326 do e = 1, n
1327 ws = 0d0
1328 !NEC$ unroll_completely
1329 do kk = 1, lx
1330 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
1331 end do
1332 w(i,j,k,e) = w(i,j,k,e) + ws
1333 end do
1334 end do
1335 end do
1336 end do
1337
1338 do j = 1, lx
1339 do i = 1, lx
1340 do k = 1, lx
1341 do e = 1, n
1342 wt = 0d0
1343 !NEC$ unroll_completely
1344 do kk = 1, lx
1345 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
1346 end do
1347 w(i,j,k,e) = w(i,j,k,e) + wt
1348 end do
1349 end do
1350 end do
1351 end do
1352
1353
1354 end subroutine sx_ax_helm_compr_lx7
1355
1356 subroutine sx_ax_helm_compr_lx6(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1357 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
1358 integer, parameter :: lx = 6
1359 integer, intent(in) :: n, m
1360 integer, intent(in) :: ind(n)
1361 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1362 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1363 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1364 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
1365 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
1366 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
1367 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
1368 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
1369 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
1370 real(kind=rp), intent(in) :: dx(lx, lx)
1371 real(kind=rp), intent(in) :: dy(lx, lx)
1372 real(kind=rp), intent(in) :: dz(lx, lx)
1373 real(kind=rp), intent(in) :: dxt(lx, lx)
1374 real(kind=rp), intent(in) :: dyt(lx, lx)
1375 real(kind=rp), intent(in) :: dzt(lx, lx)
1376 integer :: e, c, i, j, k, jj, kk
1377 real(kind=rp) :: ur(lx, lx, lx, n)
1378 real(kind=rp) :: us(lx, lx, lx, n)
1379 real(kind=rp) :: ut(lx, lx, lx, n)
1380 real(kind=rp) :: uur(lx, lx, lx, n)
1381 real(kind=rp) :: uus(lx, lx, lx, n)
1382 real(kind=rp) :: uut(lx, lx, lx, n)
1383 real(kind=rp) :: wr, ws, wt
1384
1385 do i = 1, lx
1386 do jj = 1, lx * lx * n
1387 wr = 0d0
1388 do kk = 1, lx
1389 wr = wr + dx(i,kk)*u(kk,jj,1,1)
1390 end do
1391 ur(i,jj,1,1) = wr
1392 end do
1393 end do
1394
1395 do k = 1, lx
1396 do i = 1, lx
1397 do j = 1, lx
1398 do e = 1, n
1399 ws = 0d0
1400 !NEC$ unroll_completely
1401 do kk = 1, lx
1402 ws = ws + dy(j,kk)*u(i,kk,k,e)
1403 end do
1404 us(i,j,k,e) = ws
1405 end do
1406 end do
1407 end do
1408 end do
1409
1410 do j = 1, lx
1411 do i = 1, lx
1412 do k = 1, lx
1413 do e = 1, n
1414 wt = 0d0
1415 !NEC$ unroll_completely
1416 do kk = 1, lx
1417 wt = wt + dz(k,kk)*u(i,j,kk,e)
1418 end do
1419 ut(i,j,k,e) = wt
1420 end do
1421 end do
1422 end do
1423 end do
1424
1425 do i = 1, lx * lx * lx
1426 do e = 1, n
1427 c = ind(e)
1428 uur(i,1,1,e) = h1(i,1,1,e) * &
1429 ( g11(i,1,1,c) * ur(i,1,1,e) &
1430 + g12(i,1,1,c) * us(i,1,1,e) &
1431 + g13(i,1,1,c) * ut(i,1,1,e))
1432
1433 uus(i,1,1,e) = h1(i,1,1,e) * &
1434 ( g22(i,1,1,c) * us(i,1,1,e) &
1435 + g12(i,1,1,c) * ur(i,1,1,e) &
1436 + g23(i,1,1,c) * ut(i,1,1,e) )
1437
1438 uut(i,1,1,e) = h1(i,1,1,e) * &
1439 ( g33(i,1,1,c) * ut(i,1,1,e) &
1440 + g13(i,1,1,c) * ur(i,1,1,e) &
1441 + g23(i,1,1,c) * us(i,1,1,e))
1442 end do
1443 end do
1444
1445 do i = 1, lx
1446 do jj = 1, lx * lx * n
1447 wr = 0d0
1448 do kk = 1, lx
1449 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
1450 end do
1451 w(i,jj,1,1) = wr
1452 end do
1453 end do
1454
1455 do k = 1, lx
1456 do i = 1, lx
1457 do j = 1, lx
1458 do e = 1, n
1459 ws = 0d0
1460 !NEC$ unroll_completely
1461 do kk = 1, lx
1462 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
1463 end do
1464 w(i,j,k,e) = w(i,j,k,e) + ws
1465 end do
1466 end do
1467 end do
1468 end do
1469
1470 do j = 1, lx
1471 do i = 1, lx
1472 do k = 1, lx
1473 do e = 1, n
1474 wt = 0d0
1475 !NEC$ unroll_completely
1476 do kk = 1, lx
1477 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
1478 end do
1479 w(i,j,k,e) = w(i,j,k,e) + wt
1480 end do
1481 end do
1482 end do
1483 end do
1484
1485
1486 end subroutine sx_ax_helm_compr_lx6
1487
1488 subroutine sx_ax_helm_compr_lx5(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1489 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
1490 integer, parameter :: lx = 5
1491 integer, intent(in) :: n, m
1492 integer, intent(in) :: ind(n)
1493 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1494 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1495 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1496 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
1497 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
1498 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
1499 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
1500 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
1501 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
1502 real(kind=rp), intent(in) :: dx(lx, lx)
1503 real(kind=rp), intent(in) :: dy(lx, lx)
1504 real(kind=rp), intent(in) :: dz(lx, lx)
1505 real(kind=rp), intent(in) :: dxt(lx, lx)
1506 real(kind=rp), intent(in) :: dyt(lx, lx)
1507 real(kind=rp), intent(in) :: dzt(lx, lx)
1508 integer :: e, c, i, j, k, jj, kk
1509 real(kind=rp) :: ur(lx, lx, lx, n)
1510 real(kind=rp) :: us(lx, lx, lx, n)
1511 real(kind=rp) :: ut(lx, lx, lx, n)
1512 real(kind=rp) :: uur(lx, lx, lx, n)
1513 real(kind=rp) :: uus(lx, lx, lx, n)
1514 real(kind=rp) :: uut(lx, lx, lx, n)
1515 real(kind=rp) :: wr, ws, wt
1516
1517 do i = 1, lx
1518 do jj = 1, lx * lx * n
1519 wr = 0d0
1520 do kk = 1, lx
1521 wr = wr + dx(i,kk)*u(kk,jj,1,1)
1522 end do
1523 ur(i,jj,1,1) = wr
1524 end do
1525 end do
1526
1527 do k = 1, lx
1528 do i = 1, lx
1529 do j = 1, lx
1530 do e = 1, n
1531 ws = 0d0
1532 !NEC$ unroll_completely
1533 do kk = 1, lx
1534 ws = ws + dy(j,kk)*u(i,kk,k,e)
1535 end do
1536 us(i,j,k,e) = ws
1537 end do
1538 end do
1539 end do
1540 end do
1541
1542 do j = 1, lx
1543 do i = 1, lx
1544 do k = 1, lx
1545 do e = 1, n
1546 wt = 0d0
1547 !NEC$ unroll_completely
1548 do kk = 1, lx
1549 wt = wt + dz(k,kk)*u(i,j,kk,e)
1550 end do
1551 ut(i,j,k,e) = wt
1552 end do
1553 end do
1554 end do
1555 end do
1556
1557 do i = 1, lx * lx * lx
1558 do e = 1, n
1559 c = ind(e)
1560 uur(i,1,1,e) = h1(i,1,1,e) * &
1561 ( g11(i,1,1,c) * ur(i,1,1,e) &
1562 + g12(i,1,1,c) * us(i,1,1,e) &
1563 + g13(i,1,1,c) * ut(i,1,1,e))
1564
1565 uus(i,1,1,e) = h1(i,1,1,e) * &
1566 ( g22(i,1,1,c) * us(i,1,1,e) &
1567 + g12(i,1,1,c) * ur(i,1,1,e) &
1568 + g23(i,1,1,c) * ut(i,1,1,e) )
1569
1570 uut(i,1,1,e) = h1(i,1,1,e) * &
1571 ( g33(i,1,1,c) * ut(i,1,1,e) &
1572 + g13(i,1,1,c) * ur(i,1,1,e) &
1573 + g23(i,1,1,c) * us(i,1,1,e))
1574 end do
1575 end do
1576
1577 do i = 1, lx
1578 do jj = 1, lx * lx * n
1579 wr = 0d0
1580 do kk = 1, lx
1581 wr = wr + dxt(i,kk) * uur(kk,jj,1,1)
1582 end do
1583 w(i,jj,1,1) = wr
1584 end do
1585 end do
1586
1587 do k = 1, lx
1588 do i = 1, lx
1589 do j = 1, lx
1590 do e = 1, n
1591 ws = 0d0
1592 !NEC$ unroll_completely
1593 do kk = 1, lx
1594 ws = ws + dyt(j, kk)*uus(i,kk,k,e)
1595 end do
1596 w(i,j,k,e) = w(i,j,k,e) + ws
1597 end do
1598 end do
1599 end do
1600 end do
1601
1602 do j = 1, lx
1603 do i = 1, lx
1604 do k = 1, lx
1605 do e = 1, n
1606 wt = 0d0
1607 !NEC$ unroll_completely
1608 do kk = 1, lx
1609 wt = wt + dzt(k, kk)*uut(i,j,kk,e)
1610 end do
1611 w(i,j,k,e) = w(i,j,k,e) + wt
1612 end do
1613 end do
1614 end do
1615 end do
1616
1617
1618 end subroutine sx_ax_helm_compr_lx5
1619
1620 subroutine sx_ax_helm_compr_lx4(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1621 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
1622 integer, parameter :: lx = 4
1623 integer, intent(in) :: n, m
1624 integer, intent(in) :: ind(n)
1625 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1626 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1627 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1628 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
1629 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
1630 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
1631 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
1632 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
1633 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
1634 real(kind=rp), intent(in) :: dx(lx, lx)
1635 real(kind=rp), intent(in) :: dy(lx, lx)
1636 real(kind=rp), intent(in) :: dz(lx, lx)
1637 real(kind=rp), intent(in) :: dxt(lx, lx)
1638 real(kind=rp), intent(in) :: dyt(lx, lx)
1639 real(kind=rp), intent(in) :: dzt(lx, lx)
1640 integer :: e, c, i, j, k, jj
1641 real(kind=rp) :: ur(lx, lx, lx, n)
1642 real(kind=rp) :: us(lx, lx, lx, n)
1643 real(kind=rp) :: ut(lx, lx, lx, n)
1644 real(kind=rp) :: uur(lx, lx, lx, n)
1645 real(kind=rp) :: uus(lx, lx, lx, n)
1646 real(kind=rp) :: uut(lx, lx, lx, n)
1647
1648 do i = 1, lx
1649 do jj = 1, lx * lx * n
1650 ur(i,jj,1,1) = dx(i,1)*u(1,jj,1,1) &
1651 + dx(i,2)*u(2,jj,1,1) &
1652 + dx(i,3)*u(3,jj,1,1) &
1653 + dx(i,4)*u(4,jj,1,1)
1654 end do
1655 end do
1656
1657
1658 do k = 1, lx
1659 do i = 1, lx
1660 do j = 1, lx
1661 do e = 1, n
1662 us(i,j,k,e) = dy(j,1) * u(i,1,k,e) &
1663 + dy(j,2) * u(i,2,k,e) &
1664 + dy(j,3) * u(i,3,k,e) &
1665 + dy(j,4) * u(i,4,k,e)
1666 end do
1667 end do
1668 end do
1669 end do
1670
1671
1672
1673 do j = 1, lx
1674 do i = 1, lx
1675 do k = 1, lx
1676 do e = 1, n
1677 ut(i,j,k,e) = dz(k,1) * u(i,j,1,e) &
1678 + dz(k,2) * u(i,j,2,e) &
1679 + dz(k,3) * u(i,j,3,e) &
1680 + dz(k,4) * u(i,j,4,e)
1681 end do
1682 end do
1683 end do
1684 end do
1685
1686 do i = 1, lx * lx * lx
1687 do e = 1, n
1688 c = ind(e)
1689 uur(i,1,1,e) = h1(i,1,1,e) * &
1690 ( g11(i,1,1,c) * ur(i,1,1,e) &
1691 + g12(i,1,1,c) * us(i,1,1,e) &
1692 + g13(i,1,1,c) * ut(i,1,1,e))
1693
1694 uus(i,1,1,e) = h1(i,1,1,e) * &
1695 ( g22(i,1,1,c) * us(i,1,1,e) &
1696 + g12(i,1,1,c) * ur(i,1,1,e) &
1697 + g23(i,1,1,c) * ut(i,1,1,e) )
1698
1699 uut(i,1,1,e) = h1(i,1,1,e) * &
1700 ( g33(i,1,1,c) * ut(i,1,1,e) &
1701 + g13(i,1,1,c) * ur(i,1,1,e) &
1702 + g23(i,1,1,c) * us(i,1,1,e))
1703 end do
1704 end do
1705
1706 do i = 1, lx
1707 do jj = 1, lx * lx * n
1708 w(i,jj,1,1) = dxt(i,1) * uur(1,jj,1,1) &
1709 + dxt(i,2) * uur(2,jj,1,1) &
1710 + dxt(i,3) * uur(3,jj,1,1) &
1711 + dxt(i,4) * uur(4,jj,1,1)
1712 end do
1713 end do
1714
1715 do k = 1, lx
1716 do i = 1, lx
1717 do j = 1, lx
1718 do e = 1, n
1719 w(i,j,k,e) = w(i,j,k,e) + dyt(j,1) * uus(i,1,k,e) &
1720 + dyt(j,2) * uus(i,2,k,e) &
1721 + dyt(j,3) * uus(i,3,k,e) &
1722 + dyt(j,4) * uus(i,4,k,e)
1723 end do
1724 end do
1725 end do
1726 end do
1727
1728 do j = 1, lx
1729 do i = 1, lx
1730 do k = 1, lx
1731 do e = 1, n
1732 w(i,j,k,e) = w(i,j,k,e) + dzt(k,1) * uut(i,j,1,e) &
1733 + dzt(k,2) * uut(i,j,2,e) &
1734 + dzt(k,3) * uut(i,j,3,e) &
1735 + dzt(k,4) * uut(i,j,4,e)
1736 end do
1737 end do
1738 end do
1739 end do
1740
1741 end subroutine sx_ax_helm_compr_lx4
1742
1743 subroutine sx_ax_helm_compr_lx3(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1744 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
1745 integer, parameter :: lx = 3
1746 integer, intent(in) :: n, m
1747 integer, intent(in) :: ind(n)
1748 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1749 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1750 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1751 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
1752 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
1753 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
1754 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
1755 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
1756 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
1757 real(kind=rp), intent(in) :: dx(lx, lx)
1758 real(kind=rp), intent(in) :: dy(lx, lx)
1759 real(kind=rp), intent(in) :: dz(lx, lx)
1760 real(kind=rp), intent(in) :: dxt(lx, lx)
1761 real(kind=rp), intent(in) :: dyt(lx, lx)
1762 real(kind=rp), intent(in) :: dzt(lx, lx)
1763 integer :: e, c, i, j, k, jj
1764 real(kind=rp) :: ur(lx, lx, lx, n)
1765 real(kind=rp) :: us(lx, lx, lx, n)
1766 real(kind=rp) :: ut(lx, lx, lx, n)
1767 real(kind=rp) :: uur(lx, lx, lx, n)
1768 real(kind=rp) :: uus(lx, lx, lx, n)
1769 real(kind=rp) :: uut(lx, lx, lx, n)
1770
1771 do i = 1, lx
1772 do jj = 1, lx * lx * n
1773 ur(i,jj,1,1) = dx(i,1)*u(1,jj,1,1) &
1774 + dx(i,2)*u(2,jj,1,1) &
1775 + dx(i,3)*u(3,jj,1,1)
1776 end do
1777 end do
1778
1779
1780 do k = 1, lx
1781 do i = 1, lx
1782 do j = 1, lx
1783 do e = 1, n
1784 us(i,j,k,e) = dy(j,1) * u(i,1,k,e) &
1785 + dy(j,2) * u(i,2,k,e) &
1786 + dy(j,3) * u(i,3,k,e)
1787 end do
1788 end do
1789 end do
1790 end do
1791
1792
1793
1794 do j = 1, lx
1795 do i = 1, lx
1796 do k = 1, lx
1797 do e = 1, n
1798 ut(i,j,k,e) = dz(k,1) * u(i,j,1,e) &
1799 + dz(k,2) * u(i,j,2,e) &
1800 + dz(k,3) * u(i,j,3,e)
1801 end do
1802 end do
1803 end do
1804 end do
1805
1806 do i = 1, lx * lx * lx
1807 do e = 1, n
1808 c = ind(e)
1809 uur(i,1,1,e) = h1(i,1,1,e) * &
1810 ( g11(i,1,1,c) * ur(i,1,1,e) &
1811 + g12(i,1,1,c) * us(i,1,1,e) &
1812 + g13(i,1,1,c) * ut(i,1,1,e))
1813
1814 uus(i,1,1,e) = h1(i,1,1,e) * &
1815 ( g22(i,1,1,c) * us(i,1,1,e) &
1816 + g12(i,1,1,c) * ur(i,1,1,e) &
1817 + g23(i,1,1,c) * ut(i,1,1,e) )
1818
1819 uut(i,1,1,e) = h1(i,1,1,e) * &
1820 ( g33(i,1,1,c) * ut(i,1,1,e) &
1821 + g13(i,1,1,c) * ur(i,1,1,e) &
1822 + g23(i,1,1,c) * us(i,1,1,e))
1823 end do
1824 end do
1825
1826 do i = 1, lx
1827 do jj = 1, lx * lx * n
1828 w(i,jj,1,1) = dxt(i,1) * uur(1,jj,1,1) &
1829 + dxt(i,2) * uur(2,jj,1,1) &
1830 + dxt(i,3) * uur(3,jj,1,1)
1831 end do
1832 end do
1833
1834 do k = 1, lx
1835 do i = 1, lx
1836 do j = 1, lx
1837 do e = 1, n
1838 w(i,j,k,e) = w(i,j,k,e) + dyt(j,1) * uus(i,1,k,e) &
1839 + dyt(j,2) * uus(i,2,k,e) &
1840 + dyt(j,3) * uus(i,3,k,e)
1841 end do
1842 end do
1843 end do
1844 end do
1845
1846 do j = 1, lx
1847 do i = 1, lx
1848 do k = 1, lx
1849 do e = 1, n
1850 w(i,j,k,e) = w(i,j,k,e) + dzt(k,1) * uut(i,j,1,e) &
1851 + dzt(k,2) * uut(i,j,2,e) &
1852 + dzt(k,3) * uut(i,j,3,e)
1853 end do
1854 end do
1855 end do
1856 end do
1857
1858 end subroutine sx_ax_helm_compr_lx3
1859
1860 subroutine sx_ax_helm_compr_lx2(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1861 h1, G11, G22, G33, G12, G13, G23, ind, n, m)
1862 integer, parameter :: lx = 2
1863 integer, intent(in) :: n, m
1864 integer, intent(in) :: ind(n)
1865 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1866 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1867 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1868 real(kind=rp), intent(in) :: g11(lx, lx, lx, m)
1869 real(kind=rp), intent(in) :: g22(lx, lx, lx, m)
1870 real(kind=rp), intent(in) :: g33(lx, lx, lx, m)
1871 real(kind=rp), intent(in) :: g12(lx, lx, lx, m)
1872 real(kind=rp), intent(in) :: g13(lx, lx, lx, m)
1873 real(kind=rp), intent(in) :: g23(lx, lx, lx, m)
1874 real(kind=rp), intent(in) :: dx(lx, lx)
1875 real(kind=rp), intent(in) :: dy(lx, lx)
1876 real(kind=rp), intent(in) :: dz(lx, lx)
1877 real(kind=rp), intent(in) :: dxt(lx, lx)
1878 real(kind=rp), intent(in) :: dyt(lx, lx)
1879 real(kind=rp), intent(in) :: dzt(lx, lx)
1880 integer :: e, c, i, j, k, jj
1881 real(kind=rp) :: ur(lx, lx, lx, n)
1882 real(kind=rp) :: us(lx, lx, lx, n)
1883 real(kind=rp) :: ut(lx, lx, lx, n)
1884 real(kind=rp) :: uur(lx, lx, lx, n)
1885 real(kind=rp) :: uus(lx, lx, lx, n)
1886 real(kind=rp) :: uut(lx, lx, lx, n)
1887
1888 do i = 1, lx
1889 do jj = 1, lx * lx * n
1890 ur(i,jj,1,1) = dx(i,1) * u(1,jj,1,1) &
1891 + dx(i,2) * u(2,jj,1,1)
1892 end do
1893 end do
1894
1895 do k = 1, lx
1896 do i = 1, lx
1897 do j = 1, lx
1898 do e = 1, n
1899 us(i,j,k,e) = dy(j,1) * u(i,1,k,e) &
1900 + dy(j,2) * u(i,2,k,e)
1901 end do
1902 end do
1903 end do
1904 end do
1905
1906 do j = 1, lx
1907 do i = 1, lx
1908 do k = 1, lx
1909 do e = 1, n
1910 ut(i,j,k,e) = dz(k,1) * u(i,j,1,e) &
1911 + dz(k,2) * u(i,j,2,e)
1912 end do
1913 end do
1914 end do
1915 end do
1916
1917 do i = 1, lx * lx * lx
1918 do e = 1, n
1919 c = ind(e)
1920 uur(i,1,1,e) = h1(i,1,1,e) * &
1921 ( g11(i,1,1,c) * ur(i,1,1,e) &
1922 + g12(i,1,1,c) * us(i,1,1,e) &
1923 + g13(i,1,1,c) * ut(i,1,1,e))
1924
1925 uus(i,1,1,e) = h1(i,1,1,e) * &
1926 ( g22(i,1,1,c) * us(i,1,1,e) &
1927 + g12(i,1,1,c) * ur(i,1,1,e) &
1928 + g23(i,1,1,c) * ut(i,1,1,e) )
1929
1930 uut(i,1,1,e) = h1(i,1,1,e) * &
1931 ( g33(i,1,1,c) * ut(i,1,1,e) &
1932 + g13(i,1,1,c) * ur(i,1,1,e) &
1933 + g23(i,1,1,c) * us(i,1,1,e))
1934 end do
1935 end do
1936
1937 do i = 1, lx
1938 do jj = 1, lx * lx * n
1939 w(i,jj,1,1) = dxt(i,1) * uur(1,jj,1,1) &
1940 + dxt(i,2) * uur(2,jj,1,1)
1941 end do
1942 end do
1943
1944 do k = 1, lx
1945 do i = 1, lx
1946 do j = 1, lx
1947 do e = 1, n
1948 w(i,j,k,e) = w(i,j,k,e) + dyt(j,1) * uus(i,1,k,e) &
1949 + dyt(j,2) * uus(i,2,k,e)
1950 end do
1951 end do
1952 end do
1953 end do
1954
1955 do j = 1, lx
1956 do i = 1, lx
1957 do k = 1, lx
1958 do e = 1, n
1959 w(i,j,k,e) = w(i,j,k,e) + dzt(k,1) * uut(i,j,1,e) &
1960 + dzt(k,2) * uut(i,j,2,e)
1961 end do
1962 end do
1963 end do
1964 end do
1965
1966 end subroutine sx_ax_helm_compr_lx2
1967
1968end module ax_helm_compr_sx
utils::neko_error
Definition utils.f90:46
ax_helm_compr_sx
SX Helmholtz matrix-vector product on compressed geometric factors.
Definition ax_helm_compr_sx.f90:36
ax_helm_compr_sx::sx_ax_helm_compr_lx3
subroutine sx_ax_helm_compr_lx3(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:1745
ax_helm_compr_sx::sx_ax_helm_compr_lx9
subroutine sx_ax_helm_compr_lx9(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:962
ax_helm_compr_sx::sx_ax_helm_compr_lx4
subroutine sx_ax_helm_compr_lx4(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:1622
ax_helm_compr_sx::sx_ax_helm_compr_lx8
subroutine sx_ax_helm_compr_lx8(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:1094
ax_helm_compr_sx::sx_ax_helm_compr_lx7
subroutine sx_ax_helm_compr_lx7(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:1226
ax_helm_compr_sx::sx_ax_helm_compr_lx11
subroutine sx_ax_helm_compr_lx11(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:698
ax_helm_compr_sx::sx_ax_helm_compr_lx
subroutine sx_ax_helm_compr_lx(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m, lx)
Definition ax_helm_compr_sx.f90:171
ax_helm_compr_sx::sx_ax_helm_compr_lx10
subroutine sx_ax_helm_compr_lx10(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:830
ax_helm_compr_sx::ax_helm_compr_sx_compute
subroutine ax_helm_compr_sx_compute(this, w, u, coef, msh, xh)
Definition ax_helm_compr_sx.f90:60
ax_helm_compr_sx::sx_ax_helm_compr_lx2
subroutine sx_ax_helm_compr_lx2(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:1862
ax_helm_compr_sx::sx_ax_helm_compr_lx6
subroutine sx_ax_helm_compr_lx6(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:1358
ax_helm_compr_sx::sx_ax_helm_compr_lx12
subroutine sx_ax_helm_compr_lx12(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:566
ax_helm_compr_sx::sx_ax_helm_compr_lx5
subroutine sx_ax_helm_compr_lx5(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:1490
ax_helm_compr_sx::sx_ax_helm_compr_lx13
subroutine sx_ax_helm_compr_lx13(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:434
ax_helm_compr_sx::sx_ax_helm_compr_lx14
subroutine sx_ax_helm_compr_lx14(w, u, dx, dy, dz, dxt, dyt, dzt, h1, g11, g22, g33, g12, g13, g23, ind, n, m)
Definition ax_helm_compr_sx.f90:302
ax_helm_sx
Definition ax_helm_sx.f90:33
coefs
Coefficients.
Definition coef.f90:34
math
Definition math.f90:60
math::addcol4
subroutine, public addcol4(a, b, c, d, n)
Returns .
Definition math.f90:1219
mesh
Defines a mesh.
Definition mesh.f90:34
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:14
space
Defines a function space.
Definition space.f90:34
utils
Utilities.
Definition utils.f90:35
ax_helm_compr_sx::ax_helm_compr_sx_t
SX matrix-vector product for a Helmholtz problem, reading the geometric factors from their compresse...
Definition ax_helm_compr_sx.f90:52
ax_helm_sx::ax_helm_sx_t
Definition ax_helm_sx.f90:43
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:136
mesh::mesh_t
Definition mesh.f90:73
space::space_t
The function space for the SEM solution fields.
Definition space.f90:64