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