Neko 1.99.3
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
pc_jacobi.f90
Go to the documentation of this file.
1! Copyright (c) 2020-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!
34module jacobi
35 use math, only : col2, col3, invcol1, addcol3
36 use precon, only : pc_t
37 use coefs, only : coef_t
38 use num_types, only : rp
39 use dofmap, only : dofmap_t
40 use gather_scatter, only : gs_t, gs_op_add
41 implicit none
42 private
43
45 type, public, extends(pc_t) :: jacobi_t
46 real(kind=rp), allocatable :: d(:, :, :, :)
47 type(gs_t), pointer :: gs_h
48 type(dofmap_t), pointer :: dof
49 type(coef_t), pointer :: coef
50 contains
51 procedure, pass(this) :: init => jacobi_init
52 procedure, pass(this) :: free => jacobi_free
53 procedure, pass(this) :: solve => jacobi_solve
54 procedure, pass(this) :: update => jacobi_update
55 end type jacobi_t
56
57contains
58
59 subroutine jacobi_init(this, coef, dof, gs_h)
60 class(jacobi_t), intent(inout) :: this
61 type(coef_t), intent(in), target :: coef
62 type(dofmap_t), intent(in), target :: dof
63 type(gs_t), intent(inout), target :: gs_h
64
65 call this%free()
66 this%gs_h => gs_h
67 this%dof => dof
68 this%coef => coef
69 allocate(this%d(dof%Xh%lx, dof%Xh%ly, dof%Xh%lz, dof%msh%nelv))
70 call jacobi_update(this)
71
72 end subroutine jacobi_init
73
74 subroutine jacobi_free(this)
75 class(jacobi_t), intent(inout) :: this
76 if (allocated(this%d)) then
77 deallocate(this%d)
78 end if
79 nullify(this%dof)
80 nullify(this%gs_h)
81 nullify(this%coef)
82 end subroutine jacobi_free
83
86 subroutine jacobi_solve(this, z, r, n)
87 integer, intent(in) :: n
88 class(jacobi_t), intent(inout) :: this
89 real(kind=rp), dimension(n), intent(inout) :: z
90 real(kind=rp), dimension(n), intent(inout) :: r
91 call col3(z, r, this%d, n)
92 end subroutine jacobi_solve
93
95 subroutine jacobi_update(this)
96 class(jacobi_t), intent(inout) :: this
97 associate(dof => this%dof, coef => this%coef, gs_h => this%gs_h)
98
99 !$omp parallel
100 select case (dof%Xh%lx)
101 case (14)
102 call jacobi_update_lx14(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
103 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
104 dof%msh%dfrmd_el, dof%msh%nelv)
105 case (13)
106 call jacobi_update_lx13(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
107 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
108 dof%msh%dfrmd_el, dof%msh%nelv)
109 case (12)
110 call jacobi_update_lx12(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
111 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
112 dof%msh%dfrmd_el, dof%msh%nelv)
113 case (11)
114 call jacobi_update_lx11(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
115 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
116 dof%msh%dfrmd_el, dof%msh%nelv)
117 case (10)
118 call jacobi_update_lx10(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
119 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
120 dof%msh%dfrmd_el, dof%msh%nelv)
121 case (9)
122 call jacobi_update_lx9(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
123 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
124 dof%msh%dfrmd_el, dof%msh%nelv)
125 case (8)
126 call jacobi_update_lx8(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
127 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
128 dof%msh%dfrmd_el, dof%msh%nelv)
129 case (7)
130 call jacobi_update_lx7(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
131 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
132 dof%msh%dfrmd_el, dof%msh%nelv)
133 case (6)
134 call jacobi_update_lx6(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
135 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
136 dof%msh%dfrmd_el, dof%msh%nelv)
137 case (5)
138 call jacobi_update_lx5(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
139 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
140 dof%msh%dfrmd_el, dof%msh%nelv)
141 case (4)
142 call jacobi_update_lx4(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
143 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
144 dof%msh%dfrmd_el, dof%msh%nelv)
145 case (3)
146 call jacobi_update_lx3(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
147 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
148 dof%msh%dfrmd_el, dof%msh%nelv)
149 case (2)
150 call jacobi_update_lx2(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
151 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
152 dof%msh%dfrmd_el, dof%msh%nelv)
153 case default
154 call jacobi_update_lx(this%d, dof%Xh%dxt, dof%Xh%dyt, dof%Xh%dzt, &
155 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
156 dof%msh%dfrmd_el, dof%msh%nelv, dof%Xh%lx)
157 end select
158 !$omp end parallel
159
160 call col2(this%d, coef%h1, coef%dof%size())
161 if (coef%ifh2) call addcol3(this%d, coef%h2, coef%B, coef%dof%size())
162 call gs_h%op(this%d, dof%size(), gs_op_add)
163 call invcol1(this%d, dof%size())
164 end associate
165 end subroutine jacobi_update
166
168 subroutine jacobi_update_lx(d, dxt, dyt, dzt, G11, G22, G33, &
169 G12, G13, G23, dfrmd_el, n, lx)
170 integer, intent(in) :: n, lx
171 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
172 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
173 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
174 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
175 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
176 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
177 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
178 real(kind=rp), intent(in) :: dxt(lx, lx)
179 real(kind=rp), intent(in) :: dyt(lx, lx)
180 real(kind=rp), intent(in) :: dzt(lx, lx)
181 logical, intent(in) :: dfrmd_el(n)
182 integer :: i, j, k, l, e
183
184 !$omp do
185 do e = 1, n
186 do k = 1, lx
187 do j = 1, lx
188 !$omp simd
189 do i = 1, lx*lx*lx
190 d(i,j,k,e) = 0.0_rp
191 end do
192 end do
193 end do
194
195 do l = 1, lx
196 do k = 1, lx
197 do j = 1, lx
198 !$omp simd
199 do i = 1, lx
200 d(i,j,k,e) = d(i,j,k,e) + &
201 g11(l,j,k,e) * dxt(i,l)**2
202 end do
203 end do
204 end do
205 end do
206 do l = 1, lx
207 do k = 1, lx
208 do j = 1, lx
209 !$omp simd
210 do i = 1, lx
211 d(i,j,k,e) = d(i,j,k,e) + &
212 g22(i,l,k,e) * dyt(j,l)**2
213 end do
214 end do
215 end do
216 end do
217 do l = 1, lx
218 do k = 1, lx
219 do j = 1, lx
220 !$omp simd
221 do i = 1, lx
222 d(i,j,k,e) = d(i,j,k,e) + &
223 g33(i,j,l,e) * dzt(k,l)**2
224 end do
225 end do
226 end do
227 end do
228
229 if (dfrmd_el(e)) then
230 do j = 1, lx, lx-1
231 do k = 1, lx, lx-1
232 d(1,j,k,e) = d(1,j,k,e) &
233 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
234 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
235 d(lx,j,k,e) = d(lx,j,k,e) &
236 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
237 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
238 end do
239 end do
240
241 do i = 1, lx, lx-1
242 do k = 1, lx, lx-1
243 d(i,1,k,e) = d(i,1,k,e) &
244 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
245 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
246 d(i,lx,k,e) = d(i,lx,k,e) &
247 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
248 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
249 end do
250 end do
251 do i = 1, lx, lx-1
252 do j = 1, lx, lx-1
253 d(i,j,1,e) = d(i,j,1,e) &
254 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
255 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
256 d(i,j,lx,e) = d(i,j,lx,e) &
257 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
258 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
259 end do
260 end do
261 end if
262 end do
263 !$omp end do
264 end subroutine jacobi_update_lx
265
266 subroutine jacobi_update_lx14(d, dxt, dyt, dzt, G11, G22, G33, &
267 G12, G13, G23, dfrmd_el, n)
268 integer, parameter :: lx = 14
269 integer, intent(in) :: n
270 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
271 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
272 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
273 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
274 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
275 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
276 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
277 real(kind=rp), intent(in) :: dxt(lx, lx)
278 real(kind=rp), intent(in) :: dyt(lx, lx)
279 real(kind=rp), intent(in) :: dzt(lx, lx)
280 logical, intent(in) :: dfrmd_el(n)
281 integer :: i, j, k, l, e
282
283 !$omp do
284 do e = 1, n
285 do k = 1, lx
286 do j = 1, lx
287 !$omp simd
288 do i = 1, lx
289 d(i,j,k,e) = 0.0_rp
290 end do
291 end do
292 end do
293
294 do l = 1, lx
295 do k = 1, lx
296 do j = 1, lx
297 !$omp simd
298 do i = 1, lx
299 d(i,j,k,e) = d(i,j,k,e) + &
300 g11(l,j,k,e) * dxt(i,l)**2
301 end do
302 end do
303 end do
304 end do
305 do l = 1, lx
306 do k = 1, lx
307 do j = 1, lx
308 !$omp simd
309 do i = 1, lx
310 d(i,j,k,e) = d(i,j,k,e) + &
311 g22(i,l,k,e) * dyt(j,l)**2
312 end do
313 end do
314 end do
315 end do
316 do l = 1, lx
317 do k = 1, lx
318 do j = 1, lx
319 !$omp simd
320 do i = 1, lx
321 d(i,j,k,e) = d(i,j,k,e) + &
322 g33(i,j,l,e) * dzt(k,l)**2
323 end do
324 end do
325 end do
326 end do
327
328 if (dfrmd_el(e)) then
329 do j = 1, lx, lx-1
330 do k = 1, lx, lx-1
331 d(1,j,k,e) = d(1,j,k,e) &
332 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
333 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
334 d(lx,j,k,e) = d(lx,j,k,e) &
335 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
336 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
337 end do
338 end do
339
340 do i = 1, lx, lx-1
341 do k = 1, lx, lx-1
342 d(i,1,k,e) = d(i,1,k,e) &
343 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
344 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
345 d(i,lx,k,e) = d(i,lx,k,e) &
346 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
347 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
348 end do
349 end do
350 do i = 1, lx, lx-1
351 do j = 1, lx, lx-1
352 d(i,j,1,e) = d(i,j,1,e) &
353 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
354 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
355 d(i,j,lx,e) = d(i,j,lx,e) &
356 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
357 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
358 end do
359 end do
360 end if
361 end do
362 !$omp end do
363 end subroutine jacobi_update_lx14
364
365 subroutine jacobi_update_lx13(d, dxt, dyt, dzt, G11, G22, G33, &
366 G12, G13, G23, dfrmd_el, n)
367 integer, parameter :: lx = 13
368 integer, intent(in) :: n
369 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
370 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
371 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
372 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
373 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
374 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
375 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
376 real(kind=rp), intent(in) :: dxt(lx, lx)
377 real(kind=rp), intent(in) :: dyt(lx, lx)
378 real(kind=rp), intent(in) :: dzt(lx, lx)
379 logical, intent(in) :: dfrmd_el(n)
380 integer :: i, j, k, l, e
381
382 !$omp do
383 do e = 1, n
384 do k = 1, lx
385 do j = 1, lx
386 !$omp simd
387 do i = 1, lx
388 d(i,j,k,e) = 0.0_rp
389 end do
390 end do
391 end do
392
393 do l = 1, lx
394 do k = 1, lx
395 do j = 1, lx
396 !$omp simd
397 do i = 1, lx
398 d(i,j,k,e) = d(i,j,k,e) + &
399 g11(l,j,k,e) * dxt(i,l)**2
400 end do
401 end do
402 end do
403 end do
404 do l = 1, lx
405 do k = 1, lx
406 do j = 1, lx
407 !$omp simd
408 do i = 1, lx
409 d(i,j,k,e) = d(i,j,k,e) + &
410 g22(i,l,k,e) * dyt(j,l)**2
411 end do
412 end do
413 end do
414 end do
415 do l = 1, lx
416 do k = 1, lx
417 do j = 1, lx
418 !$omp simd
419 do i = 1, lx
420 d(i,j,k,e) = d(i,j,k,e) + &
421 g33(i,j,l,e) * dzt(k,l)**2
422 end do
423 end do
424 end do
425 end do
426
427 if (dfrmd_el(e)) then
428 do j = 1, lx, lx-1
429 do k = 1, lx, lx-1
430 d(1,j,k,e) = d(1,j,k,e) &
431 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
432 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
433 d(lx,j,k,e) = d(lx,j,k,e) &
434 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
435 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
436 end do
437 end do
438
439 do i = 1, lx, lx-1
440 do k = 1, lx, lx-1
441 d(i,1,k,e) = d(i,1,k,e) &
442 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
443 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
444 d(i,lx,k,e) = d(i,lx,k,e) &
445 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
446 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
447 end do
448 end do
449 do i = 1, lx, lx-1
450 do j = 1, lx, lx-1
451 d(i,j,1,e) = d(i,j,1,e) &
452 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
453 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
454 d(i,j,lx,e) = d(i,j,lx,e) &
455 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
456 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
457 end do
458 end do
459 end if
460 end do
461 !$omp end do
462 end subroutine jacobi_update_lx13
463
464 subroutine jacobi_update_lx12(d, dxt, dyt, dzt, G11, G22, G33, &
465 G12, G13, G23, dfrmd_el, n)
466 integer, parameter :: lx = 12
467 integer, intent(in) :: n
468 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
469 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
470 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
471 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
472 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
473 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
474 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
475 real(kind=rp), intent(in) :: dxt(lx, lx)
476 real(kind=rp), intent(in) :: dyt(lx, lx)
477 real(kind=rp), intent(in) :: dzt(lx, lx)
478 logical, intent(in) :: dfrmd_el(n)
479 integer :: i, j, k, l, e
480
481 !$omp do
482 do e = 1, n
483 do k = 1, lx
484 do j = 1, lx
485 !$omp simd
486 do i = 1, lx
487 d(i,j,k,e) = 0.0_rp
488 end do
489 end do
490 end do
491
492 do l = 1, lx
493 do k = 1, lx
494 do j = 1, lx
495 !$omp simd
496 do i = 1, lx
497 d(i,j,k,e) = d(i,j,k,e) + &
498 g11(l,j,k,e) * dxt(i,l)**2
499 end do
500 end do
501 end do
502 end do
503 do l = 1, lx
504 do k = 1, lx
505 do j = 1, lx
506 !$omp simd
507 do i = 1, lx
508 d(i,j,k,e) = d(i,j,k,e) + &
509 g22(i,l,k,e) * dyt(j,l)**2
510 end do
511 end do
512 end do
513 end do
514 do l = 1, lx
515 do k = 1, lx
516 do j = 1, lx
517 !$omp simd
518 do i = 1, lx
519 d(i,j,k,e) = d(i,j,k,e) + &
520 g33(i,j,l,e) * dzt(k,l)**2
521 end do
522 end do
523 end do
524 end do
525
526 if (dfrmd_el(e)) then
527 do j = 1, lx, lx-1
528 do k = 1, lx, lx-1
529 d(1,j,k,e) = d(1,j,k,e) &
530 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
531 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
532 d(lx,j,k,e) = d(lx,j,k,e) &
533 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
534 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
535 end do
536 end do
537
538 do i = 1, lx, lx-1
539 do k = 1, lx, lx-1
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,lx,k,e) = d(i,lx,k,e) &
544 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
545 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
546 end do
547 end do
548 do i = 1, lx, lx-1
549 do j = 1, lx, lx-1
550 d(i,j,1,e) = d(i,j,1,e) &
551 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
552 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
553 d(i,j,lx,e) = d(i,j,lx,e) &
554 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
555 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
556 end do
557 end do
558 end if
559 end do
560 !$omp end do
561 end subroutine jacobi_update_lx12
562
563 subroutine jacobi_update_lx11(d, dxt, dyt, dzt, G11, G22, G33, &
564 G12, G13, G23, dfrmd_el, n)
565 integer, parameter :: lx = 11
566 integer, intent(in) :: n
567 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
568 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
569 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
570 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
571 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
572 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
573 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
574 real(kind=rp), intent(in) :: dxt(lx, lx)
575 real(kind=rp), intent(in) :: dyt(lx, lx)
576 real(kind=rp), intent(in) :: dzt(lx, lx)
577 logical, intent(in) :: dfrmd_el(n)
578 integer :: i, j, k, l, e
579
580 !$omp do
581 do e = 1, n
582 do k = 1, lx
583 do j = 1, lx
584 !$omp simd
585 do i = 1, lx
586 d(i,j,k,e) = 0.0_rp
587 end do
588 end do
589 end do
590
591 do l = 1, lx
592 do k = 1, lx
593 do j = 1, lx
594 !$omp simd
595 do i = 1, lx
596 d(i,j,k,e) = d(i,j,k,e) + &
597 g11(l,j,k,e) * dxt(i,l)**2
598 end do
599 end do
600 end do
601 end do
602 do l = 1, lx
603 do k = 1, lx
604 do j = 1, lx
605 !$omp simd
606 do i = 1, lx
607 d(i,j,k,e) = d(i,j,k,e) + &
608 g22(i,l,k,e) * dyt(j,l)**2
609 end do
610 end do
611 end do
612 end do
613 do l = 1, lx
614 do k = 1, lx
615 do j = 1, lx
616 !$omp simd
617 do i = 1, lx
618 d(i,j,k,e) = d(i,j,k,e) + &
619 g33(i,j,l,e) * dzt(k,l)**2
620 end do
621 end do
622 end do
623 end do
624
625 if (dfrmd_el(e)) then
626 do j = 1, lx, lx-1
627 do k = 1, lx, lx-1
628 d(1,j,k,e) = d(1,j,k,e) &
629 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
630 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
631 d(lx,j,k,e) = d(lx,j,k,e) &
632 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
633 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
634 end do
635 end do
636
637 do i = 1, lx, lx-1
638 do k = 1, lx, lx-1
639 d(i,1,k,e) = d(i,1,k,e) &
640 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
641 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
642 d(i,lx,k,e) = d(i,lx,k,e) &
643 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
644 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
645 end do
646 end do
647 do i = 1, lx, lx-1
648 do j = 1, lx, lx-1
649 d(i,j,1,e) = d(i,j,1,e) &
650 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
651 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
652 d(i,j,lx,e) = d(i,j,lx,e) &
653 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
654 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
655 end do
656 end do
657 end if
658 end do
659 !$omp end do
660 end subroutine jacobi_update_lx11
661
662 subroutine jacobi_update_lx10(d, dxt, dyt, dzt, G11, G22, G33, &
663 G12, G13, G23, dfrmd_el, n)
664 integer, parameter :: lx = 10
665 integer, intent(in) :: n
666 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
667 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
668 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
669 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
670 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
671 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
672 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
673 real(kind=rp), intent(in) :: dxt(lx, lx)
674 real(kind=rp), intent(in) :: dyt(lx, lx)
675 real(kind=rp), intent(in) :: dzt(lx, lx)
676 logical, intent(in) :: dfrmd_el(n)
677 integer :: i, j, k, l, e
678
679 !$omp do
680 do e = 1, n
681 do k = 1, lx
682 do j = 1, lx
683 !$omp simd
684 do i = 1, lx
685 d(i,j,k,e) = 0.0_rp
686 end do
687 end do
688 end do
689
690 do l = 1, lx
691 do k = 1, lx
692 do j = 1, lx
693 !$omp simd
694 do i = 1, lx
695 d(i,j,k,e) = d(i,j,k,e) + &
696 g11(l,j,k,e) * dxt(i,l)**2
697 end do
698 end do
699 end do
700 end do
701 do l = 1, lx
702 do k = 1, lx
703 do j = 1, lx
704 !$omp simd
705 do i = 1, lx
706 d(i,j,k,e) = d(i,j,k,e) + &
707 g22(i,l,k,e) * dyt(j,l)**2
708 end do
709 end do
710 end do
711 end do
712 do l = 1, lx
713 do k = 1, lx
714 do j = 1, lx
715 !$omp simd
716 do i = 1, lx
717 d(i,j,k,e) = d(i,j,k,e) + &
718 g33(i,j,l,e) * dzt(k,l)**2
719 end do
720 end do
721 end do
722 end do
723
724 if (dfrmd_el(e)) then
725 do j = 1, lx, lx-1
726 do k = 1, lx, lx-1
727 d(1,j,k,e) = d(1,j,k,e) &
728 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
729 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
730 d(lx,j,k,e) = d(lx,j,k,e) &
731 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
732 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
733 end do
734 end do
735
736 do i = 1, lx, lx-1
737 do k = 1, lx, lx-1
738 d(i,1,k,e) = d(i,1,k,e) &
739 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
740 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
741 d(i,lx,k,e) = d(i,lx,k,e) &
742 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
743 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
744 end do
745 end do
746 do i = 1, lx, lx-1
747 do j = 1, lx, lx-1
748 d(i,j,1,e) = d(i,j,1,e) &
749 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
750 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
751 d(i,j,lx,e) = d(i,j,lx,e) &
752 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
753 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
754 end do
755 end do
756 end if
757 end do
758 !$omp end do
759 end subroutine jacobi_update_lx10
760
761 subroutine jacobi_update_lx9(d, dxt, dyt, dzt, G11, G22, G33, &
762 G12, G13, G23, dfrmd_el, n)
763 integer, parameter :: lx = 9
764 integer, intent(in) :: n
765 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
766 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
767 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
768 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
769 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
770 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
771 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
772 real(kind=rp), intent(in) :: dxt(lx, lx)
773 real(kind=rp), intent(in) :: dyt(lx, lx)
774 real(kind=rp), intent(in) :: dzt(lx, lx)
775 logical, intent(in) :: dfrmd_el(n)
776 integer :: i, j, k, l, e
777
778 !$omp do
779 do e = 1, n
780 do k = 1, lx
781 do j = 1, lx
782 !$omp simd
783 do i = 1, lx
784 d(i,j,k,e) = 0.0_rp
785 end do
786 end do
787 end do
788
789 do l = 1, lx
790 do k = 1, lx
791 do j = 1, lx
792 !$omp simd
793 do i = 1, lx
794 d(i,j,k,e) = d(i,j,k,e) + &
795 g11(l,j,k,e) * dxt(i,l)**2
796 end do
797 end do
798 end do
799 end do
800 do l = 1, lx
801 do k = 1, lx
802 do j = 1, lx
803 !$omp simd
804 do i = 1, lx
805 d(i,j,k,e) = d(i,j,k,e) + &
806 g22(i,l,k,e) * dyt(j,l)**2
807 end do
808 end do
809 end do
810 end do
811 do l = 1, lx
812 do k = 1, lx
813 do j = 1, lx
814 !$omp simd
815 do i = 1, lx
816 d(i,j,k,e) = d(i,j,k,e) + &
817 g33(i,j,l,e) * dzt(k,l)**2
818 end do
819 end do
820 end do
821 end do
822
823 if (dfrmd_el(e)) then
824 do j = 1, lx, lx-1
825 do k = 1, lx, lx-1
826 d(1,j,k,e) = d(1,j,k,e) &
827 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
828 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
829 d(lx,j,k,e) = d(lx,j,k,e) &
830 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
831 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
832 end do
833 end do
834
835 do i = 1, lx, lx-1
836 do k = 1, lx, lx-1
837 d(i,1,k,e) = d(i,1,k,e) &
838 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
839 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
840 d(i,lx,k,e) = d(i,lx,k,e) &
841 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
842 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
843 end do
844 end do
845 do i = 1, lx, lx-1
846 do j = 1, lx, lx-1
847 d(i,j,1,e) = d(i,j,1,e) &
848 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
849 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
850 d(i,j,lx,e) = d(i,j,lx,e) &
851 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
852 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
853 end do
854 end do
855 end if
856 end do
857 !$omp end do
858 end subroutine jacobi_update_lx9
859
860 subroutine jacobi_update_lx8(d, dxt, dyt, dzt, G11, G22, G33, &
861 G12, G13, G23, dfrmd_el, n)
862 integer, parameter :: lx = 8
863 integer, intent(in) :: n
864 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
865 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
866 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
867 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
868 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
869 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
870 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
871 real(kind=rp), intent(in) :: dxt(lx, lx)
872 real(kind=rp), intent(in) :: dyt(lx, lx)
873 real(kind=rp), intent(in) :: dzt(lx, lx)
874 logical, intent(in) :: dfrmd_el(n)
875 integer :: i, j, k, l, e
876
877 !$omp do
878 do e = 1, n
879 do k = 1, lx
880 do j = 1, lx
881 !$omp simd
882 do i = 1, lx
883 d(i,j,k,e) = 0.0_rp
884 end do
885 end do
886 end do
887
888 do l = 1, lx
889 do k = 1, lx
890 do j = 1, lx
891 !$omp simd
892 do i = 1, lx
893 d(i,j,k,e) = d(i,j,k,e) + &
894 g11(l,j,k,e) * dxt(i,l)**2
895 end do
896 end do
897 end do
898 end do
899 do l = 1, lx
900 do k = 1, lx
901 do j = 1, lx
902 !$omp simd
903 do i = 1, lx
904 d(i,j,k,e) = d(i,j,k,e) + &
905 g22(i,l,k,e) * dyt(j,l)**2
906 end do
907 end do
908 end do
909 end do
910 do l = 1, lx
911 do k = 1, lx
912 do j = 1, lx
913 !$omp simd
914 do i = 1, lx
915 d(i,j,k,e) = d(i,j,k,e) + &
916 g33(i,j,l,e) * dzt(k,l)**2
917 end do
918 end do
919 end do
920 end do
921
922 if (dfrmd_el(e)) then
923 do j = 1, lx, lx-1
924 do k = 1, lx, lx-1
925 d(1,j,k,e) = d(1,j,k,e) &
926 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
927 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
928 d(lx,j,k,e) = d(lx,j,k,e) &
929 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
930 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
931 end do
932 end do
933
934 do i = 1, lx, lx-1
935 do k = 1, lx, lx-1
936 d(i,1,k,e) = d(i,1,k,e) &
937 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
938 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
939 d(i,lx,k,e) = d(i,lx,k,e) &
940 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
941 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
942 end do
943 end do
944 do i = 1, lx, lx-1
945 do j = 1, lx, lx-1
946 d(i,j,1,e) = d(i,j,1,e) &
947 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
948 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
949 d(i,j,lx,e) = d(i,j,lx,e) &
950 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
951 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
952 end do
953 end do
954 end if
955 end do
956 !$omp end do
957 end subroutine jacobi_update_lx8
958
959 subroutine jacobi_update_lx7(d, dxt, dyt, dzt, G11, G22, G33, &
960 G12, G13, G23, dfrmd_el, n)
961 integer, parameter :: lx = 7
962 integer, intent(in) :: n
963 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
964 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
965 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
966 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
967 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
968 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
969 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
970 real(kind=rp), intent(in) :: dxt(lx, lx)
971 real(kind=rp), intent(in) :: dyt(lx, lx)
972 real(kind=rp), intent(in) :: dzt(lx, lx)
973 logical, intent(in) :: dfrmd_el(n)
974 integer :: i, j, k, l, e
975
976 !$omp do
977 do e = 1, n
978 do k = 1, lx
979 do j = 1, lx
980 !$omp simd
981 do i = 1, lx
982 d(i,j,k,e) = 0.0_rp
983 end do
984 end do
985 end do
986
987 do l = 1, lx
988 do k = 1, lx
989 do j = 1, lx
990 !$omp simd
991 do i = 1, lx
992 d(i,j,k,e) = d(i,j,k,e) + &
993 g11(l,j,k,e) * dxt(i,l)**2
994 end do
995 end do
996 end do
997 end do
998 do l = 1, lx
999 do k = 1, lx
1000 do j = 1, lx
1001 !$omp simd
1002 do i = 1, lx
1003 d(i,j,k,e) = d(i,j,k,e) + &
1004 g22(i,l,k,e) * dyt(j,l)**2
1005 end do
1006 end do
1007 end do
1008 end do
1009 do l = 1, lx
1010 do k = 1, lx
1011 do j = 1, lx
1012 !$omp simd
1013 do i = 1, lx
1014 d(i,j,k,e) = d(i,j,k,e) + &
1015 g33(i,j,l,e) * dzt(k,l)**2
1016 end do
1017 end do
1018 end do
1019 end do
1020
1021 if (dfrmd_el(e)) then
1022 do j = 1, lx, lx-1
1023 do k = 1, lx, lx-1
1024 d(1,j,k,e) = d(1,j,k,e) &
1025 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1026 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1027 d(lx,j,k,e) = d(lx,j,k,e) &
1028 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1029 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1030 end do
1031 end do
1032
1033 do i = 1, lx, lx-1
1034 do k = 1, lx, lx-1
1035 d(i,1,k,e) = d(i,1,k,e) &
1036 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1037 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1038 d(i,lx,k,e) = d(i,lx,k,e) &
1039 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1040 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1041 end do
1042 end do
1043 do i = 1, lx, lx-1
1044 do j = 1, lx, lx-1
1045 d(i,j,1,e) = d(i,j,1,e) &
1046 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1047 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1048 d(i,j,lx,e) = d(i,j,lx,e) &
1049 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1050 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1051 end do
1052 end do
1053 end if
1054 end do
1055 !$omp end do
1056 end subroutine jacobi_update_lx7
1057
1058 subroutine jacobi_update_lx6(d, dxt, dyt, dzt, G11, G22, G33, &
1059 G12, G13, G23, dfrmd_el, n)
1060 integer, parameter :: lx = 6
1061 integer, intent(in) :: n
1062 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1063 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1064 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1065 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1066 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1067 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1068 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1069 real(kind=rp), intent(in) :: dxt(lx, lx)
1070 real(kind=rp), intent(in) :: dyt(lx, lx)
1071 real(kind=rp), intent(in) :: dzt(lx, lx)
1072 logical, intent(in) :: dfrmd_el(n)
1073 integer :: i, j, k, l, e
1074
1075 !$omp do
1076 do e = 1, n
1077 do k = 1, lx
1078 do j = 1, lx
1079 !$omp simd
1080 do i = 1, lx
1081 d(i,j,k,e) = 0.0_rp
1082 end do
1083 end do
1084 end do
1085
1086 do l = 1, lx
1087 do k = 1, lx
1088 do j = 1, lx
1089 !$omp simd
1090 do i = 1, lx
1091 d(i,j,k,e) = d(i,j,k,e) + &
1092 g11(l,j,k,e) * dxt(i,l)**2
1093 end do
1094 end do
1095 end do
1096 end do
1097 do l = 1, lx
1098 do k = 1, lx
1099 do j = 1, lx
1100 !$omp simd
1101 do i = 1, lx
1102 d(i,j,k,e) = d(i,j,k,e) + &
1103 g22(i,l,k,e) * dyt(j,l)**2
1104 end do
1105 end do
1106 end do
1107 end do
1108 do l = 1, lx
1109 do k = 1, lx
1110 do j = 1, lx
1111 !$omp simd
1112 do i = 1, lx
1113 d(i,j,k,e) = d(i,j,k,e) + &
1114 g33(i,j,l,e) * dzt(k,l)**2
1115 end do
1116 end do
1117 end do
1118 end do
1119
1120 if (dfrmd_el(e)) then
1121 do j = 1, lx, lx-1
1122 do k = 1, lx, lx-1
1123 d(1,j,k,e) = d(1,j,k,e) &
1124 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1125 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1126 d(lx,j,k,e) = d(lx,j,k,e) &
1127 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1128 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1129 end do
1130 end do
1131
1132 do i = 1, lx, lx-1
1133 do k = 1, lx, lx-1
1134 d(i,1,k,e) = d(i,1,k,e) &
1135 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1136 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1137 d(i,lx,k,e) = d(i,lx,k,e) &
1138 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1139 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1140 end do
1141 end do
1142 do i = 1, lx, lx-1
1143 do j = 1, lx, lx-1
1144 d(i,j,1,e) = d(i,j,1,e) &
1145 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1146 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1147 d(i,j,lx,e) = d(i,j,lx,e) &
1148 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1149 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1150 end do
1151 end do
1152 end if
1153 end do
1154 !$omp end do
1155 end subroutine jacobi_update_lx6
1156
1157 subroutine jacobi_update_lx5(d, dxt, dyt, dzt, G11, G22, G33, &
1158 G12, G13, G23, dfrmd_el, n)
1159 integer, parameter :: lx = 5
1160 integer, intent(in) :: n
1161 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1162 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1163 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1164 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1165 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1166 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1167 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1168 real(kind=rp), intent(in) :: dxt(lx, lx)
1169 real(kind=rp), intent(in) :: dyt(lx, lx)
1170 real(kind=rp), intent(in) :: dzt(lx, lx)
1171 logical, intent(in) :: dfrmd_el(n)
1172 integer :: i, j, k, l, e
1173
1174 !$omp do
1175 do e = 1, n
1176 do k = 1, lx
1177 do j = 1, lx
1178 !$omp simd
1179 do i = 1, lx
1180 d(i,j,k,e) = 0.0_rp
1181 end do
1182 end do
1183 end do
1184
1185 do l = 1, lx
1186 do k = 1, lx
1187 do j = 1, lx
1188 !$omp simd
1189 do i = 1, lx
1190 d(i,j,k,e) = d(i,j,k,e) + &
1191 g11(l,j,k,e) * dxt(i,l)**2
1192 end do
1193 end do
1194 end do
1195 end do
1196 do l = 1, lx
1197 do k = 1, lx
1198 do j = 1, lx
1199 !$omp simd
1200 do i = 1, lx
1201 d(i,j,k,e) = d(i,j,k,e) + &
1202 g22(i,l,k,e) * dyt(j,l)**2
1203 end do
1204 end do
1205 end do
1206 end do
1207 do l = 1, lx
1208 do k = 1, lx
1209 do j = 1, lx
1210 !$omp simd
1211 do i = 1, lx
1212 d(i,j,k,e) = d(i,j,k,e) + &
1213 g33(i,j,l,e) * dzt(k,l)**2
1214 end do
1215 end do
1216 end do
1217 end do
1218
1219 if (dfrmd_el(e)) then
1220 do j = 1, lx, lx-1
1221 do k = 1, lx, lx-1
1222 d(1,j,k,e) = d(1,j,k,e) &
1223 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1224 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1225 d(lx,j,k,e) = d(lx,j,k,e) &
1226 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1227 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1228 end do
1229 end do
1230
1231 do i = 1, lx, lx-1
1232 do k = 1, lx, lx-1
1233 d(i,1,k,e) = d(i,1,k,e) &
1234 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1235 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1236 d(i,lx,k,e) = d(i,lx,k,e) &
1237 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1238 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1239 end do
1240 end do
1241 do i = 1, lx, lx-1
1242 do j = 1, lx, lx-1
1243 d(i,j,1,e) = d(i,j,1,e) &
1244 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1245 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1246 d(i,j,lx,e) = d(i,j,lx,e) &
1247 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1248 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1249 end do
1250 end do
1251 end if
1252 end do
1253 !$omp end do
1254 end subroutine jacobi_update_lx5
1255
1256 subroutine jacobi_update_lx4(d, dxt, dyt, dzt, G11, G22, G33, &
1257 G12, G13, G23, dfrmd_el, n)
1258 integer, parameter :: lx = 4
1259 integer, intent(in) :: n
1260 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1261 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1262 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1263 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1264 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1265 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1266 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1267 real(kind=rp), intent(in) :: dxt(lx, lx)
1268 real(kind=rp), intent(in) :: dyt(lx, lx)
1269 real(kind=rp), intent(in) :: dzt(lx, lx)
1270 logical, intent(in) :: dfrmd_el(n)
1271 integer :: i, j, k, l, e
1272
1273 !$omp do
1274 do e = 1, n
1275 do k = 1, lx
1276 do j = 1, lx
1277 !$omp simd
1278 do i = 1, lx
1279 d(i,j,k,e) = 0.0_rp
1280 end do
1281 end do
1282 end do
1283
1284 do l = 1, lx
1285 do k = 1, lx
1286 do j = 1, lx
1287 !$omp simd
1288 do i = 1, lx
1289 d(i,j,k,e) = d(i,j,k,e) + &
1290 g11(l,j,k,e) * dxt(i,l)**2
1291 end do
1292 end do
1293 end do
1294 end do
1295 do l = 1, lx
1296 do k = 1, lx
1297 do j = 1, lx
1298 !$omp simd
1299 do i = 1, lx
1300 d(i,j,k,e) = d(i,j,k,e) + &
1301 g22(i,l,k,e) * dyt(j,l)**2
1302 end do
1303 end do
1304 end do
1305 end do
1306 do l = 1, lx
1307 do k = 1, lx
1308 do j = 1, lx
1309 !$omp simd
1310 do i = 1, lx
1311 d(i,j,k,e) = d(i,j,k,e) + &
1312 g33(i,j,l,e) * dzt(k,l)**2
1313 end do
1314 end do
1315 end do
1316 end do
1317
1318 if (dfrmd_el(e)) then
1319 do j = 1, lx, lx-1
1320 do k = 1, lx, lx-1
1321 d(1,j,k,e) = d(1,j,k,e) &
1322 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1323 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1324 d(lx,j,k,e) = d(lx,j,k,e) &
1325 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1326 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1327 end do
1328 end do
1329
1330 do i = 1, lx, lx-1
1331 do k = 1, lx, lx-1
1332 d(i,1,k,e) = d(i,1,k,e) &
1333 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1334 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1335 d(i,lx,k,e) = d(i,lx,k,e) &
1336 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1337 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1338 end do
1339 end do
1340 do i = 1, lx, lx-1
1341 do j = 1, lx, lx-1
1342 d(i,j,1,e) = d(i,j,1,e) &
1343 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1344 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1345 d(i,j,lx,e) = d(i,j,lx,e) &
1346 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1347 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1348 end do
1349 end do
1350 end if
1351 end do
1352 !$omp end do
1353 end subroutine jacobi_update_lx4
1354
1355 subroutine jacobi_update_lx3(d, dxt, dyt, dzt, G11, G22, G33, &
1356 G12, G13, G23, dfrmd_el, n)
1357 integer, parameter :: lx = 3
1358 integer, intent(in) :: n
1359 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1360 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1361 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1362 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1363 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1364 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1365 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1366 real(kind=rp), intent(in) :: dxt(lx, lx)
1367 real(kind=rp), intent(in) :: dyt(lx, lx)
1368 real(kind=rp), intent(in) :: dzt(lx, lx)
1369 logical, intent(in) :: dfrmd_el(n)
1370 integer :: i, j, k, l, e
1371
1372 !$omp do
1373 do e = 1, n
1374 do k = 1, lx
1375 do j = 1, lx
1376 !$omp simd
1377 do i = 1, lx
1378 d(i,j,k,e) = 0.0_rp
1379 end do
1380 end do
1381 end do
1382
1383 do l = 1, lx
1384 do k = 1, lx
1385 do j = 1, lx
1386 !$omp simd
1387 do i = 1, lx
1388 d(i,j,k,e) = d(i,j,k,e) + &
1389 g11(l,j,k,e) * dxt(i,l)**2
1390 end do
1391 end do
1392 end do
1393 end do
1394 do l = 1, lx
1395 do k = 1, lx
1396 do j = 1, lx
1397 !$omp simd
1398 do i = 1, lx
1399 d(i,j,k,e) = d(i,j,k,e) + &
1400 g22(i,l,k,e) * dyt(j,l)**2
1401 end do
1402 end do
1403 end do
1404 end do
1405 do l = 1, lx
1406 do k = 1, lx
1407 do j = 1, lx
1408 !$omp simd
1409 do i = 1, lx
1410 d(i,j,k,e) = d(i,j,k,e) + &
1411 g33(i,j,l,e) * dzt(k,l)**2
1412 end do
1413 end do
1414 end do
1415 end do
1416
1417 if (dfrmd_el(e)) then
1418 do j = 1, lx, lx-1
1419 do k = 1, lx, lx-1
1420 d(1,j,k,e) = d(1,j,k,e) &
1421 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1422 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1423 d(lx,j,k,e) = d(lx,j,k,e) &
1424 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1425 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1426 end do
1427 end do
1428
1429 do i = 1, lx, lx-1
1430 do k = 1, lx, lx-1
1431 d(i,1,k,e) = d(i,1,k,e) &
1432 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1433 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1434 d(i,lx,k,e) = d(i,lx,k,e) &
1435 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1436 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1437 end do
1438 end do
1439 do i = 1, lx, lx-1
1440 do j = 1, lx, lx-1
1441 d(i,j,1,e) = d(i,j,1,e) &
1442 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1443 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1444 d(i,j,lx,e) = d(i,j,lx,e) &
1445 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1446 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1447 end do
1448 end do
1449 end if
1450 end do
1451 !$omp end do
1452 end subroutine jacobi_update_lx3
1453
1454 subroutine jacobi_update_lx2(d, dxt, dyt, dzt, G11, G22, G33, &
1455 G12, G13, G23, dfrmd_el, n)
1456 integer, parameter :: lx = 2
1457 integer, intent(in) :: n
1458 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1459 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1460 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1461 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1462 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1463 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1464 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1465 real(kind=rp), intent(in) :: dxt(lx, lx)
1466 real(kind=rp), intent(in) :: dyt(lx, lx)
1467 real(kind=rp), intent(in) :: dzt(lx, lx)
1468 logical, intent(in) :: dfrmd_el(n)
1469 integer :: i, j, k, l, e
1470
1471 !$omp do
1472 do e = 1, n
1473 do k = 1, lx
1474 do j = 1, lx
1475 !$omp simd
1476 do i = 1, lx
1477 d(i,j,k,e) = 0.0_rp
1478 end do
1479 end do
1480 end do
1481
1482 do l = 1, lx
1483 do k = 1, lx
1484 do j = 1, lx
1485 !$omp simd
1486 do i = 1, lx
1487 d(i,j,k,e) = d(i,j,k,e) + &
1488 g11(l,j,k,e) * dxt(i,l)**2
1489 end do
1490 end do
1491 end do
1492 end do
1493 do l = 1, lx
1494 do k = 1, lx
1495 do j = 1, lx
1496 !$omp simd
1497 do i = 1, lx
1498 d(i,j,k,e) = d(i,j,k,e) + &
1499 g22(i,l,k,e) * dyt(j,l)**2
1500 end do
1501 end do
1502 end do
1503 end do
1504 do l = 1, lx
1505 do k = 1, lx
1506 do j = 1, lx
1507 !$omp simd
1508 do i = 1, lx
1509 d(i,j,k,e) = d(i,j,k,e) + &
1510 g33(i,j,l,e) * dzt(k,l)**2
1511 end do
1512 end do
1513 end do
1514 end do
1515
1516 if (dfrmd_el(e)) then
1517 do j = 1, lx, lx-1
1518 do k = 1, lx, lx-1
1519 d(1,j,k,e) = d(1,j,k,e) &
1520 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1521 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1522 d(lx,j,k,e) = d(lx,j,k,e) &
1523 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1524 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1525 end do
1526 end do
1527
1528 do i = 1, lx, lx-1
1529 do k = 1, lx, lx-1
1530 d(i,1,k,e) = d(i,1,k,e) &
1531 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1532 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1533 d(i,lx,k,e) = d(i,lx,k,e) &
1534 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1535 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1536 end do
1537 end do
1538 do i = 1, lx, lx-1
1539 do j = 1, lx, lx-1
1540 d(i,j,1,e) = d(i,j,1,e) &
1541 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1542 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1543 d(i,j,lx,e) = d(i,j,lx,e) &
1544 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1545 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1546 end do
1547 end do
1548 end if
1549 end do
1550 !$omp end do
1551 end subroutine jacobi_update_lx2
1552
1553end module jacobi
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:126
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
jacobi
Jacobi preconditioner.
Definition pc_jacobi.f90:34
jacobi::jacobi_update_lx8
subroutine jacobi_update_lx8(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:862
jacobi::jacobi_solve
subroutine jacobi_solve(this, z, r, n)
The jacobi preconditioner where .
Definition pc_jacobi.f90:87
jacobi::jacobi_update
subroutine jacobi_update(this)
Update Jacobi preconditioner if the geometry G has changed.
Definition pc_jacobi.f90:96
jacobi::jacobi_update_lx7
subroutine jacobi_update_lx7(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:961
jacobi::jacobi_update_lx11
subroutine jacobi_update_lx11(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:565
jacobi::jacobi_update_lx10
subroutine jacobi_update_lx10(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:664
jacobi::jacobi_init
subroutine jacobi_init(this, coef, dof, gs_h)
Definition pc_jacobi.f90:60
jacobi::jacobi_free
subroutine jacobi_free(this)
Definition pc_jacobi.f90:75
jacobi::jacobi_update_lx
subroutine jacobi_update_lx(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n, lx)
Generic CPU kernel for updating the Jacobi preconditioner.
Definition pc_jacobi.f90:170
jacobi::jacobi_update_lx13
subroutine jacobi_update_lx13(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:367
jacobi::jacobi_update_lx12
subroutine jacobi_update_lx12(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:466
jacobi::jacobi_update_lx5
subroutine jacobi_update_lx5(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:1159
jacobi::jacobi_update_lx9
subroutine jacobi_update_lx9(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:763
jacobi::jacobi_update_lx14
subroutine jacobi_update_lx14(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:268
jacobi::jacobi_update_lx4
subroutine jacobi_update_lx4(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:1258
jacobi::jacobi_update_lx3
subroutine jacobi_update_lx3(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:1357
jacobi::jacobi_update_lx6
subroutine jacobi_update_lx6(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:1060
jacobi::jacobi_update_lx2
subroutine jacobi_update_lx2(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, dfrmd_el, n)
Definition pc_jacobi.f90:1456
math
Definition math.f90:60
math::addcol3
subroutine, public addcol3(a, b, c, n)
Returns .
Definition math.f90:1162
math::invcol1
subroutine, public invcol1(a, n)
Invert a vector .
Definition math.f90:769
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:1044
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:1059
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
precon
Krylov preconditioner.
Definition precon.f90:34
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:63
dofmap::dofmap_t
Definition dofmap.f90:56
gather_scatter::gs_t
Definition gather_scatter.f90:70
jacobi::jacobi_t
Defines a jacobi preconditioner.
Definition pc_jacobi.f90:45
precon::pc_t
Defines a canonical Krylov preconditioner.
Definition precon.f90:40