Neko 1.99.9
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
pc_jacobi_sx.f90
Go to the documentation of this file.
1! Copyright (c) 2021, 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!
34module sx_jacobi
35 use math, only : col2, col3, addcol3, invcol1
36 use precon, only : pc_t
37 use bc_list, only : bc_list_t
38 use json_module, only : json_file
39 use coefs, only : coef_t
40 use dofmap, only : dofmap_t
41 use num_types, only : rp
42 use gather_scatter, only : gs_t, gs_op_add
43 implicit none
44 private
45
47 type, public, extends(pc_t) :: sx_jacobi_t
48 real(kind=rp), allocatable :: d(:, :, :, :)
49 type(gs_t), pointer :: gs_h
50 type(dofmap_t), pointer :: dof
51 type(coef_t), pointer :: coef
52 contains
53 procedure, pass(this) :: init => sx_jacobi_init
54 procedure, pass(this) :: init_from_components => &
55 sx_jacobi_init_from_components
56 procedure, pass(this) :: free => sx_jacobi_free
57 procedure, pass(this) :: solve => sx_jacobi_solve
58 procedure, pass(this) :: update => sx_jacobi_update
59 end type sx_jacobi_t
60
61contains
62
65 subroutine sx_jacobi_init(this, coef, bclst, json)
66 class(sx_jacobi_t), intent(inout), target :: this
67 type(coef_t), intent(in), target :: coef
68 type(bc_list_t), intent(inout), target :: bclst
69 type(json_file), intent(inout) :: json
70
71 call this%init_from_components(coef, coef%dof, coef%gs_h)
72
73 end subroutine sx_jacobi_init
74
75 subroutine sx_jacobi_init_from_components(this, coef, dof, gs_h)
76 class(sx_jacobi_t), intent(inout) :: this
77 type(coef_t), intent(in), target :: coef
78 type(dofmap_t), intent(in), target :: dof
79 type(gs_t), intent(in), target :: gs_h
80
81 call this%free()
82 this%gs_h => gs_h
83 this%dof => dof
84 this%coef => coef
85 allocate(this%d(dof%Xh%lx, dof%Xh%ly, dof%Xh%lz, dof%msh%nelv))
86 call sx_jacobi_update(this)
87
88 end subroutine sx_jacobi_init_from_components
89
90 subroutine sx_jacobi_free(this)
91 class(sx_jacobi_t), intent(inout) :: this
92 if (allocated(this%d)) then
93 deallocate(this%d)
94 end if
95 nullify(this%dof)
96 nullify(this%gs_h)
97 nullify(this%coef)
98 end subroutine sx_jacobi_free
99
102 subroutine sx_jacobi_solve(this, z, r, n)
103 integer, intent(in) :: n
104 class(sx_jacobi_t), intent(inout) :: this
105 real(kind=rp), dimension(n), intent(inout) :: z
106 real(kind=rp), dimension(n), intent(inout) :: r
107 call col3(z, r, this%d, n)
108 end subroutine sx_jacobi_solve
109
110 subroutine sx_jacobi_update(this)
111 class(sx_jacobi_t), intent(inout) :: this
112 integer :: lz, ly, lx
113 associate(dof => this%dof, coef => this%coef, &
114 gs_h => this%gs_h, nelv => this%dof%msh%nelv)
115
116 lx = dof%Xh%lx
117 ly = dof%Xh%ly
118 lz = dof%Xh%lz
119
120 this%d = 0d0
121
122 select case(lx)
123 case (14)
124 call sx_update_lx14(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
125 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
126 case (13)
127 call sx_update_lx13(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
128 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
129 case (12)
130 call sx_update_lx12(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
131 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
132 case (11)
133 call sx_update_lx11(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
134 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
135 case (10)
136 call sx_update_lx10(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
137 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
138 case (9)
139 call sx_update_lx9(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
140 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
141 case (8)
142 call sx_update_lx8(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
143 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
144 case (7)
145 call sx_update_lx7(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
146 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
147 case (6)
148 call sx_update_lx6(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
149 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
150 case (5)
151 call sx_update_lx5(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
152 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
153 case (4)
154 call sx_update_lx4(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
155 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
156 case (2)
157 call sx_update_lx2(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
158 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
159 case default
160 call sx_update_lx(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
161 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
162 nelv, lx)
163 end select
164
165 call col2(this%d, coef%h1, coef%dof%size())
166 if (coef%ifh2) call addcol3(this%d, coef%h2, coef%B, coef%dof%size())
167 call gs_h%op(this%d, dof%size(), gs_op_add)
168 if (.not. coef%ifh2) call col2(this%d, coef%mult, coef%dof%size())
169 call invcol1(this%d, dof%size())
170 end associate
171 end subroutine sx_jacobi_update
172
173 subroutine sx_update_lx(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n, lx)
174 integer, intent(in) :: n, lx
175 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
176 real(kind=rp), intent(in) :: dxt(lx, lx)
177 real(kind=rp), intent(in) :: dyt(lx, lx)
178 real(kind=rp), intent(in) :: dzt(lx, lx)
179 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
180 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
181 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
182 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
183 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
184 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
185 integer :: i, j, k, l, e
186
187
188 do l = 1, lx
189 do k = 1, lx
190 do j = 1, lx
191 do i = 1, lx
192 do e = 1, n
193 d(i, j, k, e) = d(i, j, k, e) + &
194 g11(l, j, k, e) * dxt(i, l)**2
195
196 d(i, j, k, e) = d(i, j, k, e) + &
197 g22(i, l, k, e) * dyt(j, l)**2
198
199 d(i, j, k, e) = d(i, j, k, e) + &
200 g33(i, j, l, e) * dzt(k, l)**2
201
202 end do
203 end do
204 end do
205 end do
206 end do
207
208 do j = 1, lx, lx-1
209 do k = 1, lx, lx-1
210 do e = 1, n
211 d(1, j, k, e) = d(1, j, k, e) &
212 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
213 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
214 d(lx, j, k, e) = d(lx, j, k, e) &
215 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
216 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
217 end do
218 end do
219 end do
220
221 do i = 1, lx, lx-1
222 do k = 1, lx, lx-1
223 do e = 1, n
224 d(i, 1, k, e) = d(i, 1, k, e) &
225 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
226 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
227 d(i, lx, k, e) = d(i, lx, k, e) &
228 + g12(i, lx, k, e) * dyt(lx, lx)*dxt(i, i) &
229 + g23(i, lx, k, e) * dyt(lx, lx)*dzt(k, k)
230 end do
231 end do
232 end do
233
234 do i = 1, lx, lx-1
235 do j = 1, lx, lx-1
236 do e = 1, n
237 d(i, j, 1, e) = d(i, j, 1, e) &
238 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
239 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
240 d(i, j, lx, e) = d(i, j, lx, e) &
241 + g13(i, j, lx, e) * dzt(lx, lx)*dxt(i, i) &
242 + g23(i, j, lx, e) * dzt(lx, lx)*dyt(j, j)
243 end do
244 end do
245 end do
246
247 end subroutine sx_update_lx
248
249 subroutine sx_update_lx14(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
250 integer, parameter :: lx = 14
251 integer, parameter :: ly = 14
252 integer, parameter :: lz = 14
253 integer, intent(in) :: n
254 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
255 real(kind=rp), intent(in) :: dxt(lx, lx)
256 real(kind=rp), intent(in) :: dyt(lx, lx)
257 real(kind=rp), intent(in) :: dzt(lx, lx)
258 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
259 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
260 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
261 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
262 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
263 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
264 integer :: i, j, k, l, e
265
266
267 do l = 1, lx
268 do k = 1, lz
269 do j = 1, ly
270 do i = 1, lx
271 do e = 1, n
272 d(i, j, k, e) = d(i, j, k, e) + &
273 g11(l, j, k, e) * dxt(i, l)**2
274
275 d(i, j, k, e) = d(i, j, k, e) + &
276 g22(i, l, k, e) * dyt(j, l)**2
277
278 d(i, j, k, e) = d(i, j, k, e) + &
279 g33(i, j, l, e) * dzt(k, l)**2
280
281 end do
282 end do
283 end do
284 end do
285 end do
286
287 do j = 1, ly, ly-1
288 do k = 1, lz, lz-1
289 do e = 1, n
290 d(1, j, k, e) = d(1, j, k, e) &
291 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
292 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
293 d(lx, j, k, e) = d(lx, j, k, e) &
294 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
295 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
296 end do
297 end do
298 end do
299
300 do i = 1, lx, lx-1
301 do k = 1, lz, lz-1
302 do e = 1, n
303 d(i, 1, k, e) = d(i, 1, k, e) &
304 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
305 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
306 d(i, ly, k, e) = d(i, ly, k, e) &
307 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
308 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
309 end do
310 end do
311 end do
312
313 do i = 1, lx, lx-1
314 do j = 1, ly, ly-1
315 do e = 1, n
316 d(i, j, 1, e) = d(i, j, 1, e) &
317 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
318 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
319 d(i, j, lz, e) = d(i, j, lz, e) &
320 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
321 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
322 end do
323 end do
324 end do
325
326 end subroutine sx_update_lx14
327
328 subroutine sx_update_lx13(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
329 integer, parameter :: lx = 13
330 integer, parameter :: ly = 13
331 integer, parameter :: lz = 13
332 integer, intent(in) :: n
333 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
334 real(kind=rp), intent(in) :: dxt(lx, lx)
335 real(kind=rp), intent(in) :: dyt(lx, lx)
336 real(kind=rp), intent(in) :: dzt(lx, lx)
337 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
338 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
339 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
340 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
341 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
342 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
343 integer :: i, j, k, l, e
344
345
346 do l = 1, lx
347 do k = 1, lz
348 do j = 1, ly
349 do i = 1, lx
350 do e = 1, n
351 d(i, j, k, e) = d(i, j, k, e) + &
352 g11(l, j, k, e) * dxt(i, l)**2
353
354 d(i, j, k, e) = d(i, j, k, e) + &
355 g22(i, l, k, e) * dyt(j, l)**2
356
357 d(i, j, k, e) = d(i, j, k, e) + &
358 g33(i, j, l, e) * dzt(k, l)**2
359
360 end do
361 end do
362 end do
363 end do
364 end do
365
366 do j = 1, ly, ly-1
367 do k = 1, lz, lz-1
368 do e = 1, n
369 d(1, j, k, e) = d(1, j, k, e) &
370 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
371 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
372 d(lx, j, k, e) = d(lx, j, k, e) &
373 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
374 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
375 end do
376 end do
377 end do
378
379 do i = 1, lx, lx-1
380 do k = 1, lz, lz-1
381 do e = 1, n
382 d(i, 1, k, e) = d(i, 1, k, e) &
383 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
384 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
385 d(i, ly, k, e) = d(i, ly, k, e) &
386 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
387 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
388 end do
389 end do
390 end do
391
392 do i = 1, lx, lx-1
393 do j = 1, ly, ly-1
394 do e = 1, n
395 d(i, j, 1, e) = d(i, j, 1, e) &
396 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
397 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
398 d(i, j, lz, e) = d(i, j, lz, e) &
399 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
400 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
401 end do
402 end do
403 end do
404
405 end subroutine sx_update_lx13
406
407 subroutine sx_update_lx12(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
408 integer, parameter :: lx = 12
409 integer, parameter :: ly = 12
410 integer, parameter :: lz = 12
411 integer, intent(in) :: n
412 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
413 real(kind=rp), intent(in) :: dxt(lx, lx)
414 real(kind=rp), intent(in) :: dyt(lx, lx)
415 real(kind=rp), intent(in) :: dzt(lx, lx)
416 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
417 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
418 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
419 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
420 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
421 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
422 integer :: i, j, k, l, e
423
424
425 do l = 1, lx
426 do k = 1, lz
427 do j = 1, ly
428 do i = 1, lx
429 do e = 1, n
430 d(i, j, k, e) = d(i, j, k, e) + &
431 g11(l, j, k, e) * dxt(i, l)**2
432
433 d(i, j, k, e) = d(i, j, k, e) + &
434 g22(i, l, k, e) * dyt(j, l)**2
435
436 d(i, j, k, e) = d(i, j, k, e) + &
437 g33(i, j, l, e) * dzt(k, l)**2
438
439 end do
440 end do
441 end do
442 end do
443 end do
444
445 do j = 1, ly, ly-1
446 do k = 1, lz, lz-1
447 do e = 1, n
448 d(1, j, k, e) = d(1, j, k, e) &
449 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
450 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
451 d(lx, j, k, e) = d(lx, j, k, e) &
452 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
453 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
454 end do
455 end do
456 end do
457
458 do i = 1, lx, lx-1
459 do k = 1, lz, lz-1
460 do e = 1, n
461 d(i, 1, k, e) = d(i, 1, k, e) &
462 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
463 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
464 d(i, ly, k, e) = d(i, ly, k, e) &
465 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
466 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
467 end do
468 end do
469 end do
470
471 do i = 1, lx, lx-1
472 do j = 1, ly, ly-1
473 do e = 1, n
474 d(i, j, 1, e) = d(i, j, 1, e) &
475 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
476 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
477 d(i, j, lz, e) = d(i, j, lz, e) &
478 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
479 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
480 end do
481 end do
482 end do
483
484 end subroutine sx_update_lx12
485
486 subroutine sx_update_lx11(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
487 integer, parameter :: lx = 11
488 integer, parameter :: ly = 11
489 integer, parameter :: lz = 11
490 integer, intent(in) :: n
491 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
492 real(kind=rp), intent(in) :: dxt(lx, lx)
493 real(kind=rp), intent(in) :: dyt(lx, lx)
494 real(kind=rp), intent(in) :: dzt(lx, lx)
495 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
496 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
497 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
498 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
499 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
500 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
501 integer :: i, j, k, l, e
502
503
504 do l = 1, lx
505 do k = 1, lz
506 do j = 1, ly
507 do i = 1, lx
508 do e = 1, n
509 d(i, j, k, e) = d(i, j, k, e) + &
510 g11(l, j, k, e) * dxt(i, l)**2
511
512 d(i, j, k, e) = d(i, j, k, e) + &
513 g22(i, l, k, e) * dyt(j, l)**2
514
515 d(i, j, k, e) = d(i, j, k, e) + &
516 g33(i, j, l, e) * dzt(k, l)**2
517
518 end do
519 end do
520 end do
521 end do
522 end do
523
524 do j = 1, ly, ly-1
525 do k = 1, lz, lz-1
526 do e = 1, n
527 d(1, j, k, e) = d(1, j, k, e) &
528 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
529 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
530 d(lx, j, k, e) = d(lx, j, k, e) &
531 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
532 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
533 end do
534 end do
535 end do
536
537 do i = 1, lx, lx-1
538 do k = 1, lz, lz-1
539 do e = 1, n
540 d(i, 1, k, e) = d(i, 1, k, e) &
541 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
542 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
543 d(i, ly, k, e) = d(i, ly, k, e) &
544 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
545 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
546 end do
547 end do
548 end do
549
550 do i = 1, lx, lx-1
551 do j = 1, ly, ly-1
552 do e = 1, n
553 d(i, j, 1, e) = d(i, j, 1, e) &
554 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
555 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
556 d(i, j, lz, e) = d(i, j, lz, e) &
557 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
558 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
559 end do
560 end do
561 end do
562
563 end subroutine sx_update_lx11
564
565 subroutine sx_update_lx10(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
566 integer, parameter :: lx = 10
567 integer, parameter :: ly = 10
568 integer, parameter :: lz = 10
569 integer, intent(in) :: n
570 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
571 real(kind=rp), intent(in) :: dxt(lx, lx)
572 real(kind=rp), intent(in) :: dyt(lx, lx)
573 real(kind=rp), intent(in) :: dzt(lx, lx)
574 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
575 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
576 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
577 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
578 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
579 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
580 integer :: i, j, k, l, e
581
582
583 do l = 1, lx
584 do k = 1, lz
585 do j = 1, ly
586 do i = 1, lx
587 do e = 1, n
588 d(i, j, k, e) = d(i, j, k, e) + &
589 g11(l, j, k, e) * dxt(i, l)**2
590
591 d(i, j, k, e) = d(i, j, k, e) + &
592 g22(i, l, k, e) * dyt(j, l)**2
593
594 d(i, j, k, e) = d(i, j, k, e) + &
595 g33(i, j, l, e) * dzt(k, l)**2
596
597 end do
598 end do
599 end do
600 end do
601 end do
602
603 do j = 1, ly, ly-1
604 do k = 1, lz, lz-1
605 do e = 1, n
606 d(1, j, k, e) = d(1, j, k, e) &
607 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
608 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
609 d(lx, j, k, e) = d(lx, j, k, e) &
610 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
611 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
612 end do
613 end do
614 end do
615
616 do i = 1, lx, lx-1
617 do k = 1, lz, lz-1
618 do e = 1, n
619 d(i, 1, k, e) = d(i, 1, k, e) &
620 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
621 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
622 d(i, ly, k, e) = d(i, ly, k, e) &
623 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
624 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
625 end do
626 end do
627 end do
628
629 do i = 1, lx, lx-1
630 do j = 1, ly, ly-1
631 do e = 1, n
632 d(i, j, 1, e) = d(i, j, 1, e) &
633 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
634 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
635 d(i, j, lz, e) = d(i, j, lz, e) &
636 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
637 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
638 end do
639 end do
640 end do
641
642 end subroutine sx_update_lx10
643
644 subroutine sx_update_lx9(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
645 integer, parameter :: lx = 9
646 integer, parameter :: ly = 9
647 integer, parameter :: lz = 9
648 integer, intent(in) :: n
649 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
650 real(kind=rp), intent(in) :: dxt(lx, lx)
651 real(kind=rp), intent(in) :: dyt(lx, lx)
652 real(kind=rp), intent(in) :: dzt(lx, lx)
653 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
654 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
655 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
656 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
657 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
658 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
659 integer :: i, j, k, l, e
660
661
662 do l = 1, lx
663 do k = 1, lz
664 do j = 1, ly
665 do i = 1, lx
666 do e = 1, n
667 d(i, j, k, e) = d(i, j, k, e) + &
668 g11(l, j, k, e) * dxt(i, l)**2
669
670 d(i, j, k, e) = d(i, j, k, e) + &
671 g22(i, l, k, e) * dyt(j, l)**2
672
673 d(i, j, k, e) = d(i, j, k, e) + &
674 g33(i, j, l, e) * dzt(k, l)**2
675
676 end do
677 end do
678 end do
679 end do
680 end do
681
682 do j = 1, ly, ly-1
683 do k = 1, lz, lz-1
684 do e = 1, n
685 d(1, j, k, e) = d(1, j, k, e) &
686 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
687 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
688 d(lx, j, k, e) = d(lx, j, k, e) &
689 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
690 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
691 end do
692 end do
693 end do
694
695 do i = 1, lx, lx-1
696 do k = 1, lz, lz-1
697 do e = 1, n
698 d(i, 1, k, e) = d(i, 1, k, e) &
699 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
700 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
701 d(i, ly, k, e) = d(i, ly, k, e) &
702 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
703 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
704 end do
705 end do
706 end do
707
708 do i = 1, lx, lx-1
709 do j = 1, ly, ly-1
710 do e = 1, n
711 d(i, j, 1, e) = d(i, j, 1, e) &
712 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
713 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
714 d(i, j, lz, e) = d(i, j, lz, e) &
715 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
716 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
717 end do
718 end do
719 end do
720
721 end subroutine sx_update_lx9
722
723 subroutine sx_update_lx8(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
724 integer, parameter :: lx = 8
725 integer, parameter :: ly = 8
726 integer, parameter :: lz = 8
727 integer, intent(in) :: n
728 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
729 real(kind=rp), intent(in) :: dxt(lx, lx)
730 real(kind=rp), intent(in) :: dyt(lx, lx)
731 real(kind=rp), intent(in) :: dzt(lx, lx)
732 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
733 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
734 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
735 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
736 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
737 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
738 integer :: i, j, k, l, e
739
740
741 do l = 1, lx
742 do k = 1, lz
743 do j = 1, ly
744 do i = 1, lx
745 do e = 1, n
746 d(i, j, k, e) = d(i, j, k, e) + &
747 g11(l, j, k, e) * dxt(i, l)**2
748
749 d(i, j, k, e) = d(i, j, k, e) + &
750 g22(i, l, k, e) * dyt(j, l)**2
751
752 d(i, j, k, e) = d(i, j, k, e) + &
753 g33(i, j, l, e) * dzt(k, l)**2
754
755 end do
756 end do
757 end do
758 end do
759 end do
760
761 do j = 1, ly, ly-1
762 do k = 1, lz, lz-1
763 do e = 1, n
764 d(1, j, k, e) = d(1, j, k, e) &
765 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
766 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
767 d(lx, j, k, e) = d(lx, j, k, e) &
768 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
769 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
770 end do
771 end do
772 end do
773
774 do i = 1, lx, lx-1
775 do k = 1, lz, lz-1
776 do e = 1, n
777 d(i, 1, k, e) = d(i, 1, k, e) &
778 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
779 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
780 d(i, ly, k, e) = d(i, ly, k, e) &
781 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
782 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
783 end do
784 end do
785 end do
786
787 do i = 1, lx, lx-1
788 do j = 1, ly, ly-1
789 do e = 1, n
790 d(i, j, 1, e) = d(i, j, 1, e) &
791 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
792 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
793 d(i, j, lz, e) = d(i, j, lz, e) &
794 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
795 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
796 end do
797 end do
798 end do
799
800 end subroutine sx_update_lx8
801
802 subroutine sx_update_lx7(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
803 integer, parameter :: lx = 7
804 integer, parameter :: ly = 7
805 integer, parameter :: lz = 7
806 integer, intent(in) :: n
807 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
808 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
809 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
810 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
811 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
812 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
813 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
814 real(kind=rp), intent(in) :: dxt(lx, lx)
815 real(kind=rp), intent(in) :: dyt(lx, lx)
816 real(kind=rp), intent(in) :: dzt(lx, lx)
817 integer :: i, j, k, l, e
818
819 do l = 1, lx
820 do k = 1, lz
821 do j = 1, ly
822 do i = 1, lx
823 do e = 1, n
824 d(i, j, k, e) = d(i, j, k, e) + &
825 g11(l, j, k, e) * dxt(i, l)**2
826
827 d(i, j, k, e) = d(i, j, k, e) + &
828 g22(i, l, k, e) * dyt(j, l)**2
829
830 d(i, j, k, e) = d(i, j, k, e) + &
831 g33(i, j, l, e) * dzt(k, l)**2
832
833 end do
834 end do
835 end do
836 end do
837 end do
838
839 do j = 1, ly, ly-1
840 do k = 1, lz, lz-1
841 do e = 1, n
842 d(1, j, k, e) = d(1, j, k, e) &
843 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
844 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
845 d(lx, j, k, e) = d(lx, j, k, e) &
846 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
847 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
848 end do
849 end do
850 end do
851
852 do i = 1, lx, lx-1
853 do k = 1, lz, lz-1
854 do e = 1, n
855 d(i, 1, k, e) = d(i, 1, k, e) &
856 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
857 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
858 d(i, ly, k, e) = d(i, ly, k, e) &
859 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
860 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
861 end do
862 end do
863 end do
864
865 do i = 1, lx, lx-1
866 do j = 1, ly, ly-1
867 do e = 1, n
868 d(i, j, 1, e) = d(i, j, 1, e) &
869 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
870 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
871 d(i, j, lz, e) = d(i, j, lz, e) &
872 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
873 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
874 end do
875 end do
876 end do
877
878
879 end subroutine sx_update_lx7
880
881 subroutine sx_update_lx6(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
882 integer, parameter :: lx = 6
883 integer, parameter :: ly = 6
884 integer, parameter :: lz = 6
885 integer, intent(in) :: n
886 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
887 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
888 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
889 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
890 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
891 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
892 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
893 real(kind=rp), intent(in) :: dxt(lx, lx)
894 real(kind=rp), intent(in) :: dyt(lx, lx)
895 real(kind=rp), intent(in) :: dzt(lx, lx)
896 integer :: i, j, k, l, e
897
898 do l = 1, lx
899 do k = 1, lz
900 do j = 1, ly
901 do i = 1, lx
902 do e = 1, n
903 d(i, j, k, e) = d(i, j, k, e) + &
904 g11(l, j, k, e) * dxt(i, l)**2
905
906 d(i, j, k, e) = d(i, j, k, e) + &
907 g22(i, l, k, e) * dyt(j, l)**2
908
909 d(i, j, k, e) = d(i, j, k, e) + &
910 g33(i, j, l, e) * dzt(k, l)**2
911
912 end do
913 end do
914 end do
915 end do
916 end do
917
918 do j = 1, ly, ly-1
919 do k = 1, lz, lz-1
920 do e = 1, n
921 d(1, j, k, e) = d(1, j, k, e) &
922 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
923 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
924 d(lx, j, k, e) = d(lx, j, k, e) &
925 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
926 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
927 end do
928 end do
929 end do
930
931 do i = 1, lx, lx-1
932 do k = 1, lz, lz-1
933 do e = 1, n
934 d(i, 1, k, e) = d(i, 1, k, e) &
935 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
936 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
937 d(i, ly, k, e) = d(i, ly, k, e) &
938 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
939 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
940 end do
941 end do
942 end do
943
944 do i = 1, lx, lx-1
945 do j = 1, ly, ly-1
946 do e = 1, n
947 d(i, j, 1, e) = d(i, j, 1, e) &
948 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
949 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
950 d(i, j, lz, e) = d(i, j, lz, e) &
951 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
952 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
953 end do
954 end do
955 end do
956
957
958 end subroutine sx_update_lx6
959
960 subroutine sx_update_lx5(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
961 integer, parameter :: lx = 5
962 integer, parameter :: ly = 5
963 integer, parameter :: lz = 5
964 integer, intent(in) :: n
965 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
966 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
967 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
968 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
969 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
970 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
971 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
972 real(kind=rp), intent(in) :: dxt(lx, lx)
973 real(kind=rp), intent(in) :: dyt(lx, lx)
974 real(kind=rp), intent(in) :: dzt(lx, lx)
975 integer :: i, j, k, l, e
976
977 do l = 1, lx
978 do k = 1, lz
979 do j = 1, ly
980 do i = 1, lx
981 do e = 1, n
982 d(i, j, k, e) = d(i, j, k, e) + &
983 g11(l, j, k, e) * dxt(i, l)**2
984
985 d(i, j, k, e) = d(i, j, k, e) + &
986 g22(i, l, k, e) * dyt(j, l)**2
987
988 d(i, j, k, e) = d(i, j, k, e) + &
989 g33(i, j, l, e) * dzt(k, l)**2
990
991 end do
992 end do
993 end do
994 end do
995 end do
996
997 do j = 1, ly, ly-1
998 do k = 1, lz, lz-1
999 do e = 1, n
1000 d(1, j, k, e) = d(1, j, k, e) &
1001 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
1002 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
1003 d(lx, j, k, e) = d(lx, j, k, e) &
1004 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
1005 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
1006 end do
1007 end do
1008 end do
1009
1010 do i = 1, lx, lx-1
1011 do k = 1, lz, lz-1
1012 do e = 1, n
1013 d(i, 1, k, e) = d(i, 1, k, e) &
1014 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
1015 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
1016 d(i, ly, k, e) = d(i, ly, k, e) &
1017 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
1018 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
1019 end do
1020 end do
1021 end do
1022
1023 do i = 1, lx, lx-1
1024 do j = 1, ly, ly-1
1025 do e = 1, n
1026 d(i, j, 1, e) = d(i, j, 1, e) &
1027 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
1028 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
1029 d(i, j, lz, e) = d(i, j, lz, e) &
1030 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
1031 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
1032 end do
1033 end do
1034 end do
1035
1036
1037 end subroutine sx_update_lx5
1038
1039 subroutine sx_update_lx4(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
1040 integer, parameter :: lx = 4
1041 integer, parameter :: ly = 4
1042 integer, parameter :: lz = 4
1043 integer, intent(in) :: n
1044 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1045 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1046 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1047 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1048 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1049 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1050 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1051 real(kind=rp), intent(in) :: dxt(lx, lx)
1052 real(kind=rp), intent(in) :: dyt(lx, lx)
1053 real(kind=rp), intent(in) :: dzt(lx, lx)
1054 integer :: i, j, k, l, e
1055
1056 do l = 1, lx
1057 do k = 1, lz
1058 do j = 1, ly
1059 do i = 1, lx
1060 do e = 1, n
1061 d(i, j, k, e) = d(i, j, k, e) + &
1062 g11(l, j, k, e) * dxt(i, l)**2
1063
1064 d(i, j, k, e) = d(i, j, k, e) + &
1065 g22(i, l, k, e) * dyt(j, l)**2
1066
1067 d(i, j, k, e) = d(i, j, k, e) + &
1068 g33(i, j, l, e) * dzt(k, l)**2
1069
1070 end do
1071 end do
1072 end do
1073 end do
1074 end do
1075
1076 do j = 1, ly, ly-1
1077 do k = 1, lz, lz-1
1078 do e = 1, n
1079 d(1, j, k, e) = d(1, j, k, e) &
1080 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
1081 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
1082 d(lx, j, k, e) = d(lx, j, k, e) &
1083 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
1084 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
1085 end do
1086 end do
1087 end do
1088
1089 do i = 1, lx, lx-1
1090 do k = 1, lz, lz-1
1091 do e = 1, n
1092 d(i, 1, k, e) = d(i, 1, k, e) &
1093 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
1094 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
1095 d(i, ly, k, e) = d(i, ly, k, e) &
1096 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
1097 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
1098 end do
1099 end do
1100 end do
1101
1102 do i = 1, lx, lx-1
1103 do j = 1, ly, ly-1
1104 do e = 1, n
1105 d(i, j, 1, e) = d(i, j, 1, e) &
1106 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
1107 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
1108 d(i, j, lz, e) = d(i, j, lz, e) &
1109 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
1110 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
1111 end do
1112 end do
1113 end do
1114
1115
1116 end subroutine sx_update_lx4
1117
1118 subroutine sx_update_lx3(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
1119 integer, parameter :: lx = 3
1120 integer, parameter :: ly = 3
1121 integer, parameter :: lz = 3
1122 integer, intent(in) :: n
1123 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1124 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1125 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1126 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1127 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1128 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1129 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1130 real(kind=rp), intent(in) :: dxt(lx, lx)
1131 real(kind=rp), intent(in) :: dyt(lx, lx)
1132 real(kind=rp), intent(in) :: dzt(lx, lx)
1133 integer :: i, j, k, l, e
1134
1135
1136 do l = 1, lx
1137 do k = 1, lz
1138 do j = 1, ly
1139 do i = 1, lx
1140 do e = 1, n
1141 d(i, j, k, e) = d(i, j, k, e) + &
1142 g11(l, j, k, e) * dxt(i, l)**2
1143
1144 d(i, j, k, e) = d(i, j, k, e) + &
1145 g22(i, l, k, e) * dyt(j, l)**2
1146
1147 d(i, j, k, e) = d(i, j, k, e) + &
1148 g33(i, j, l, e) * dzt(k, l)**2
1149
1150 end do
1151 end do
1152 end do
1153 end do
1154 end do
1155
1156 do j = 1, ly, ly-1
1157 do k = 1, lz, lz-1
1158 do e = 1, n
1159 d(1, j, k, e) = d(1, j, k, e) &
1160 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
1161 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
1162 d(lx, j, k, e) = d(lx, j, k, e) &
1163 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
1164 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
1165 end do
1166 end do
1167 end do
1168
1169 do i = 1, lx, lx-1
1170 do k = 1, lz, lz-1
1171 do e = 1, n
1172 d(i, 1, k, e) = d(i, 1, k, e) &
1173 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
1174 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
1175 d(i, ly, k, e) = d(i, ly, k, e) &
1176 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
1177 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
1178 end do
1179 end do
1180 end do
1181
1182 do i = 1, lx, lx-1
1183 do j = 1, ly, ly-1
1184 do e = 1, n
1185 d(i, j, 1, e) = d(i, j, 1, e) &
1186 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
1187 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
1188 d(i, j, lz, e) = d(i, j, lz, e) &
1189 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
1190 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
1191 end do
1192 end do
1193 end do
1194
1195 end subroutine sx_update_lx3
1196
1197 subroutine sx_update_lx2(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
1198 integer, parameter :: lx = 2
1199 integer, parameter :: ly = 2
1200 integer, parameter :: lz = 2
1201 integer, intent(in) :: n
1202 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1203 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1204 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1205 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1206 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1207 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1208 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1209 real(kind=rp), intent(in) :: dxt(lx, lx)
1210 real(kind=rp), intent(in) :: dyt(lx, lx)
1211 real(kind=rp), intent(in) :: dzt(lx, lx)
1212 integer :: i, j, k, l, e
1213
1214
1215 do l = 1, lx
1216 do k = 1, lz
1217 do j = 1, ly
1218 do i = 1, lx
1219 do e = 1, n
1220 d(i, j, k, e) = d(i, j, k, e) + &
1221 g11(l, j, k, e) * dxt(i, l)**2
1222
1223 d(i, j, k, e) = d(i, j, k, e) + &
1224 g22(i, l, k, e) * dyt(j, l)**2
1225
1226 d(i, j, k, e) = d(i, j, k, e) + &
1227 g33(i, j, l, e) * dzt(k, l)**2
1228
1229 end do
1230 end do
1231 end do
1232 end do
1233 end do
1234
1235 do j = 1, ly, ly-1
1236 do k = 1, lz, lz-1
1237 do e = 1, n
1238 d(1, j, k, e) = d(1, j, k, e) &
1239 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
1240 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
1241 d(lx, j, k, e) = d(lx, j, k, e) &
1242 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
1243 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
1244 end do
1245 end do
1246 end do
1247
1248 do i = 1, lx, lx-1
1249 do k = 1, lz, lz-1
1250 do e = 1, n
1251 d(i, 1, k, e) = d(i, 1, k, e) &
1252 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
1253 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
1254 d(i, ly, k, e) = d(i, ly, k, e) &
1255 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
1256 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
1257 end do
1258 end do
1259 end do
1260
1261 do i = 1, lx, lx-1
1262 do j = 1, ly, ly-1
1263 do e = 1, n
1264 d(i, j, 1, e) = d(i, j, 1, e) &
1265 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
1266 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
1267 d(i, j, lz, e) = d(i, j, lz, e) &
1268 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
1269 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
1270 end do
1271 end do
1272 end do
1273
1274 end subroutine sx_update_lx2
1275
1276end module sx_jacobi
1277
solve
__device__ T solve(const T u, const T y, const T guess, const T nu, const T kappa, const T B)
Definition spalding_kernel.h:115
bc_list
Defines a list of bc_t.
Definition bc_list.f90:34
coefs
Coefficients.
Definition coef.f90:34
dofmap
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
gather_scatter
Gather-scatter.
Definition gather_scatter.f90:34
math
Definition math.f90:60
math::addcol3
subroutine, public addcol3(a, b, c, n)
Returns .
Definition math.f90:1203
math::invcol1
subroutine, public invcol1(a, n)
Invert a vector .
Definition math.f90:810
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:1085
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:1100
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:14
precon
Krylov preconditioner.
Definition precon.f90:34
sx_jacobi
Jacobi preconditioner SX-Aurora backend.
Definition pc_jacobi_sx.f90:34
sx_jacobi::sx_update_lx3
subroutine sx_update_lx3(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:1119
sx_jacobi::sx_jacobi_free
subroutine sx_jacobi_free(this)
Definition pc_jacobi_sx.f90:91
sx_jacobi::sx_jacobi_solve
subroutine sx_jacobi_solve(this, z, r, n)
The jacobi preconditioner where .
Definition pc_jacobi_sx.f90:103
sx_jacobi::sx_update_lx5
subroutine sx_update_lx5(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:961
sx_jacobi::sx_jacobi_init
subroutine sx_jacobi_init(this, coef, bclst, json)
Initialize from the case file. The boundary conditions and the JSON dictionary are not needed by this...
Definition pc_jacobi_sx.f90:66
sx_jacobi::sx_update_lx12
subroutine sx_update_lx12(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:408
sx_jacobi::sx_update_lx11
subroutine sx_update_lx11(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:487
sx_jacobi::sx_update_lx2
subroutine sx_update_lx2(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:1198
sx_jacobi::sx_update_lx7
subroutine sx_update_lx7(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:803
sx_jacobi::sx_jacobi_init_from_components
subroutine sx_jacobi_init_from_components(this, coef, dof, gs_h)
Definition pc_jacobi_sx.f90:76
sx_jacobi::sx_update_lx13
subroutine sx_update_lx13(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:329
sx_jacobi::sx_update_lx4
subroutine sx_update_lx4(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:1040
sx_jacobi::sx_update_lx6
subroutine sx_update_lx6(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:882
sx_jacobi::sx_update_lx
subroutine sx_update_lx(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n, lx)
Definition pc_jacobi_sx.f90:174
sx_jacobi::sx_jacobi_update
subroutine sx_jacobi_update(this)
Definition pc_jacobi_sx.f90:111
sx_jacobi::sx_update_lx9
subroutine sx_update_lx9(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:645
sx_jacobi::sx_update_lx14
subroutine sx_update_lx14(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:250
sx_jacobi::sx_update_lx8
subroutine sx_update_lx8(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:724
sx_jacobi::sx_update_lx10
subroutine sx_update_lx10(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:566
bc_list::bc_list_t
A list of allocatable `bc_t`. Follows the standard interface of lists.
Definition bc_list.f90:49
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:136
dofmap::dofmap_t
Definition dofmap.f90:58
gather_scatter::gs_t
Gather-scatter kernel.
Definition gather_scatter.f90:90
precon::pc_t
Defines a canonical Krylov preconditioner.
Definition precon.f90:44
sx_jacobi::sx_jacobi_t
Defines a jacobi preconditioner for SX-Aurora.
Definition pc_jacobi_sx.f90:47