Neko 1.99.1
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-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 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
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
159 call col2(this%d,coef%h1,coef%dof%size())
160 if (coef%ifh2) call addcol3(this%d,coef%h2,coef%B,coef%dof%size())
161 call gs_h%op(this%d, dof%size(), gs_op_add)
162 call invcol1(this%d,dof%size())
163 end associate
164 end subroutine jacobi_update
165
167 subroutine jacobi_update_lx(d, dxt, dyt, dzt, G11, G22, G33, &
168 G12, G13, G23, dfrmd_el, n, lx)
169 integer, intent(in) :: n, lx
170 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
171 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
172 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
173 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
174 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
175 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
176 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
177 real(kind=rp), intent(in) :: dxt(lx, lx)
178 real(kind=rp), intent(in) :: dyt(lx, lx)
179 real(kind=rp), intent(in) :: dzt(lx, lx)
180 logical, intent(in) :: dfrmd_el(n)
181 integer :: i, j, k, l, e
182
183 d = 0d0
184
185 do e = 1,n
186 do l = 1,lx
187 do k = 1,lx
188 do j = 1,lx
189 do i = 1,lx
190 d(i,j,k,e) = d(i,j,k,e) + &
191 g11(l,j,k,e) * dxt(i,l)**2
192 end do
193 end do
194 end do
195 end do
196 do l = 1,lx
197 do k = 1,lx
198 do j = 1,lx
199 do i = 1,lx
200 d(i,j,k,e) = d(i,j,k,e) + &
201 g22(i,l,k,e) * dyt(j,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 do i = 1,lx
210 d(i,j,k,e) = d(i,j,k,e) + &
211 g33(i,j,l,e) * dzt(k,l)**2
212 end do
213 end do
214 end do
215 end do
216
217 if (dfrmd_el(e)) then
218 do j = 1,lx,lx-1
219 do k = 1,lx,lx-1
220 d(1,j,k,e) = d(1,j,k,e) &
221 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
222 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
223 d(lx,j,k,e) = d(lx,j,k,e) &
224 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
225 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
226 end do
227 end do
228
229 do i = 1,lx,lx-1
230 do k = 1,lx,lx-1
231 d(i,1,k,e) = d(i,1,k,e) &
232 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
233 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
234 d(i,lx,k,e) = d(i,lx,k,e) &
235 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
236 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
237 end do
238 end do
239 do i = 1,lx,lx-1
240 do j = 1,lx,lx-1
241 d(i,j,1,e) = d(i,j,1,e) &
242 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
243 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
244 d(i,j,lx,e) = d(i,j,lx,e) &
245 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
246 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
247 end do
248 end do
249 end if
250 end do
251 end subroutine jacobi_update_lx
252
253 subroutine jacobi_update_lx14(d, dxt, dyt, dzt, G11, G22, G33, &
254 G12, G13, G23, dfrmd_el, n)
255 integer, parameter :: lx = 14
256 integer, intent(in) :: n
257 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
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 real(kind=rp), intent(in) :: dxt(lx, lx)
265 real(kind=rp), intent(in) :: dyt(lx, lx)
266 real(kind=rp), intent(in) :: dzt(lx, lx)
267 logical, intent(in) :: dfrmd_el(n)
268 integer :: i, j, k, l, e
269
270 d = 0d0
271
272 do e = 1,n
273 do l = 1,lx
274 do k = 1,lx
275 do j = 1,lx
276 do i = 1,lx
277 d(i,j,k,e) = d(i,j,k,e) + &
278 g11(l,j,k,e) * dxt(i,l)**2
279 end do
280 end do
281 end do
282 end do
283 do l = 1,lx
284 do k = 1,lx
285 do j = 1,lx
286 do i = 1,lx
287 d(i,j,k,e) = d(i,j,k,e) + &
288 g22(i,l,k,e) * dyt(j,l)**2
289 end do
290 end do
291 end do
292 end do
293 do l = 1,lx
294 do k = 1,lx
295 do j = 1,lx
296 do i = 1,lx
297 d(i,j,k,e) = d(i,j,k,e) + &
298 g33(i,j,l,e) * dzt(k,l)**2
299 end do
300 end do
301 end do
302 end do
303
304 if (dfrmd_el(e)) then
305 do j = 1,lx,lx-1
306 do k = 1,lx,lx-1
307 d(1,j,k,e) = d(1,j,k,e) &
308 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
309 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
310 d(lx,j,k,e) = d(lx,j,k,e) &
311 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
312 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
313 end do
314 end do
315
316 do i = 1,lx,lx-1
317 do k = 1,lx,lx-1
318 d(i,1,k,e) = d(i,1,k,e) &
319 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
320 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
321 d(i,lx,k,e) = d(i,lx,k,e) &
322 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
323 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
324 end do
325 end do
326 do i = 1,lx,lx-1
327 do j = 1,lx,lx-1
328 d(i,j,1,e) = d(i,j,1,e) &
329 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
330 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
331 d(i,j,lx,e) = d(i,j,lx,e) &
332 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
333 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
334 end do
335 end do
336 end if
337 end do
338 end subroutine jacobi_update_lx14
339
340 subroutine jacobi_update_lx13(d, dxt, dyt, dzt, G11, G22, G33, &
341 G12, G13, G23, dfrmd_el, n)
342 integer, parameter :: lx = 13
343 integer, intent(in) :: n
344 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
345 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
346 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
347 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
348 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
349 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
350 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
351 real(kind=rp), intent(in) :: dxt(lx, lx)
352 real(kind=rp), intent(in) :: dyt(lx, lx)
353 real(kind=rp), intent(in) :: dzt(lx, lx)
354 logical, intent(in) :: dfrmd_el(n)
355 integer :: i, j, k, l, e
356
357 d = 0d0
358
359 do e = 1,n
360 do l = 1,lx
361 do k = 1,lx
362 do j = 1,lx
363 do i = 1,lx
364 d(i,j,k,e) = d(i,j,k,e) + &
365 g11(l,j,k,e) * dxt(i,l)**2
366 end do
367 end do
368 end do
369 end do
370 do l = 1,lx
371 do k = 1,lx
372 do j = 1,lx
373 do i = 1,lx
374 d(i,j,k,e) = d(i,j,k,e) + &
375 g22(i,l,k,e) * dyt(j,l)**2
376 end do
377 end do
378 end do
379 end do
380 do l = 1,lx
381 do k = 1,lx
382 do j = 1,lx
383 do i = 1,lx
384 d(i,j,k,e) = d(i,j,k,e) + &
385 g33(i,j,l,e) * dzt(k,l)**2
386 end do
387 end do
388 end do
389 end do
390
391 if (dfrmd_el(e)) then
392 do j = 1,lx,lx-1
393 do k = 1,lx,lx-1
394 d(1,j,k,e) = d(1,j,k,e) &
395 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
396 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
397 d(lx,j,k,e) = d(lx,j,k,e) &
398 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
399 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
400 end do
401 end do
402
403 do i = 1,lx,lx-1
404 do k = 1,lx,lx-1
405 d(i,1,k,e) = d(i,1,k,e) &
406 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
407 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
408 d(i,lx,k,e) = d(i,lx,k,e) &
409 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
410 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
411 end do
412 end do
413 do i = 1,lx,lx-1
414 do j = 1,lx,lx-1
415 d(i,j,1,e) = d(i,j,1,e) &
416 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
417 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
418 d(i,j,lx,e) = d(i,j,lx,e) &
419 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
420 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
421 end do
422 end do
423 end if
424 end do
425 end subroutine jacobi_update_lx13
426
427 subroutine jacobi_update_lx12(d, dxt, dyt, dzt, G11, G22, G33, &
428 G12, G13, G23, dfrmd_el, n)
429 integer, parameter :: lx = 12
430 integer, intent(in) :: n
431 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
432 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
433 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
434 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
435 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
436 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
437 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
438 real(kind=rp), intent(in) :: dxt(lx, lx)
439 real(kind=rp), intent(in) :: dyt(lx, lx)
440 real(kind=rp), intent(in) :: dzt(lx, lx)
441 logical, intent(in) :: dfrmd_el(n)
442 integer :: i, j, k, l, e
443
444 d = 0d0
445
446 do e = 1,n
447 do l = 1,lx
448 do k = 1,lx
449 do j = 1,lx
450 do i = 1,lx
451 d(i,j,k,e) = d(i,j,k,e) + &
452 g11(l,j,k,e) * dxt(i,l)**2
453 end do
454 end do
455 end do
456 end do
457 do l = 1,lx
458 do k = 1,lx
459 do j = 1,lx
460 do i = 1,lx
461 d(i,j,k,e) = d(i,j,k,e) + &
462 g22(i,l,k,e) * dyt(j,l)**2
463 end do
464 end do
465 end do
466 end do
467 do l = 1,lx
468 do k = 1,lx
469 do j = 1,lx
470 do i = 1,lx
471 d(i,j,k,e) = d(i,j,k,e) + &
472 g33(i,j,l,e) * dzt(k,l)**2
473 end do
474 end do
475 end do
476 end do
477
478 if (dfrmd_el(e)) then
479 do j = 1,lx,lx-1
480 do k = 1,lx,lx-1
481 d(1,j,k,e) = d(1,j,k,e) &
482 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
483 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
484 d(lx,j,k,e) = d(lx,j,k,e) &
485 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
486 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
487 end do
488 end do
489
490 do i = 1,lx,lx-1
491 do k = 1,lx,lx-1
492 d(i,1,k,e) = d(i,1,k,e) &
493 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
494 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
495 d(i,lx,k,e) = d(i,lx,k,e) &
496 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
497 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
498 end do
499 end do
500 do i = 1,lx,lx-1
501 do j = 1,lx,lx-1
502 d(i,j,1,e) = d(i,j,1,e) &
503 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
504 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
505 d(i,j,lx,e) = d(i,j,lx,e) &
506 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
507 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
508 end do
509 end do
510 end if
511 end do
512 end subroutine jacobi_update_lx12
513
514 subroutine jacobi_update_lx11(d, dxt, dyt, dzt, G11, G22, G33, &
515 G12, G13, G23, dfrmd_el, n)
516 integer, parameter :: lx = 11
517 integer, intent(in) :: n
518 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
519 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
520 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
521 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
522 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
523 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
524 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
525 real(kind=rp), intent(in) :: dxt(lx, lx)
526 real(kind=rp), intent(in) :: dyt(lx, lx)
527 real(kind=rp), intent(in) :: dzt(lx, lx)
528 logical, intent(in) :: dfrmd_el(n)
529 integer :: i, j, k, l, e
530
531 d = 0d0
532
533 do e = 1,n
534 do l = 1,lx
535 do k = 1,lx
536 do j = 1,lx
537 do i = 1,lx
538 d(i,j,k,e) = d(i,j,k,e) + &
539 g11(l,j,k,e) * dxt(i,l)**2
540 end do
541 end do
542 end do
543 end do
544 do l = 1,lx
545 do k = 1,lx
546 do j = 1,lx
547 do i = 1,lx
548 d(i,j,k,e) = d(i,j,k,e) + &
549 g22(i,l,k,e) * dyt(j,l)**2
550 end do
551 end do
552 end do
553 end do
554 do l = 1,lx
555 do k = 1,lx
556 do j = 1,lx
557 do i = 1,lx
558 d(i,j,k,e) = d(i,j,k,e) + &
559 g33(i,j,l,e) * dzt(k,l)**2
560 end do
561 end do
562 end do
563 end do
564
565 if (dfrmd_el(e)) then
566 do j = 1,lx,lx-1
567 do k = 1,lx,lx-1
568 d(1,j,k,e) = d(1,j,k,e) &
569 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
570 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
571 d(lx,j,k,e) = d(lx,j,k,e) &
572 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
573 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
574 end do
575 end do
576
577 do i = 1,lx,lx-1
578 do k = 1,lx,lx-1
579 d(i,1,k,e) = d(i,1,k,e) &
580 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
581 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
582 d(i,lx,k,e) = d(i,lx,k,e) &
583 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
584 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
585 end do
586 end do
587 do i = 1,lx,lx-1
588 do j = 1,lx,lx-1
589 d(i,j,1,e) = d(i,j,1,e) &
590 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
591 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
592 d(i,j,lx,e) = d(i,j,lx,e) &
593 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
594 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
595 end do
596 end do
597 end if
598 end do
599 end subroutine jacobi_update_lx11
600
601 subroutine jacobi_update_lx10(d, dxt, dyt, dzt, G11, G22, G33, &
602 G12, G13, G23, dfrmd_el, n)
603 integer, parameter :: lx = 10
604 integer, intent(in) :: n
605 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
606 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
607 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
608 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
609 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
610 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
611 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
612 real(kind=rp), intent(in) :: dxt(lx, lx)
613 real(kind=rp), intent(in) :: dyt(lx, lx)
614 real(kind=rp), intent(in) :: dzt(lx, lx)
615 logical, intent(in) :: dfrmd_el(n)
616 integer :: i, j, k, l, e
617
618 d = 0d0
619
620 do e = 1,n
621 do l = 1,lx
622 do k = 1,lx
623 do j = 1,lx
624 do i = 1,lx
625 d(i,j,k,e) = d(i,j,k,e) + &
626 g11(l,j,k,e) * dxt(i,l)**2
627 end do
628 end do
629 end do
630 end do
631 do l = 1,lx
632 do k = 1,lx
633 do j = 1,lx
634 do i = 1,lx
635 d(i,j,k,e) = d(i,j,k,e) + &
636 g22(i,l,k,e) * dyt(j,l)**2
637 end do
638 end do
639 end do
640 end do
641 do l = 1,lx
642 do k = 1,lx
643 do j = 1,lx
644 do i = 1,lx
645 d(i,j,k,e) = d(i,j,k,e) + &
646 g33(i,j,l,e) * dzt(k,l)**2
647 end do
648 end do
649 end do
650 end do
651
652 if (dfrmd_el(e)) then
653 do j = 1,lx,lx-1
654 do k = 1,lx,lx-1
655 d(1,j,k,e) = d(1,j,k,e) &
656 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
657 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
658 d(lx,j,k,e) = d(lx,j,k,e) &
659 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
660 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
661 end do
662 end do
663
664 do i = 1,lx,lx-1
665 do k = 1,lx,lx-1
666 d(i,1,k,e) = d(i,1,k,e) &
667 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
668 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
669 d(i,lx,k,e) = d(i,lx,k,e) &
670 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
671 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
672 end do
673 end do
674 do i = 1,lx,lx-1
675 do j = 1,lx,lx-1
676 d(i,j,1,e) = d(i,j,1,e) &
677 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
678 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
679 d(i,j,lx,e) = d(i,j,lx,e) &
680 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
681 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
682 end do
683 end do
684 end if
685 end do
686 end subroutine jacobi_update_lx10
687
688 subroutine jacobi_update_lx9(d, dxt, dyt, dzt, G11, G22, G33, &
689 G12, G13, G23, dfrmd_el, n)
690 integer, parameter :: lx = 9
691 integer, intent(in) :: n
692 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
693 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
694 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
695 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
696 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
697 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
698 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
699 real(kind=rp), intent(in) :: dxt(lx, lx)
700 real(kind=rp), intent(in) :: dyt(lx, lx)
701 real(kind=rp), intent(in) :: dzt(lx, lx)
702 logical, intent(in) :: dfrmd_el(n)
703 integer :: i, j, k, l, e
704
705 d = 0d0
706
707 do e = 1,n
708 do l = 1,lx
709 do k = 1,lx
710 do j = 1,lx
711 do i = 1,lx
712 d(i,j,k,e) = d(i,j,k,e) + &
713 g11(l,j,k,e) * dxt(i,l)**2
714 end do
715 end do
716 end do
717 end do
718 do l = 1,lx
719 do k = 1,lx
720 do j = 1,lx
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 do i = 1,lx
732 d(i,j,k,e) = d(i,j,k,e) + &
733 g33(i,j,l,e) * dzt(k,l)**2
734 end do
735 end do
736 end do
737 end do
738
739 if (dfrmd_el(e)) then
740 do j = 1,lx,lx-1
741 do k = 1,lx,lx-1
742 d(1,j,k,e) = d(1,j,k,e) &
743 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
744 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
745 d(lx,j,k,e) = d(lx,j,k,e) &
746 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
747 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
748 end do
749 end do
750
751 do i = 1,lx,lx-1
752 do k = 1,lx,lx-1
753 d(i,1,k,e) = d(i,1,k,e) &
754 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
755 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
756 d(i,lx,k,e) = d(i,lx,k,e) &
757 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
758 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
759 end do
760 end do
761 do i = 1,lx,lx-1
762 do j = 1,lx,lx-1
763 d(i,j,1,e) = d(i,j,1,e) &
764 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
765 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
766 d(i,j,lx,e) = d(i,j,lx,e) &
767 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
768 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
769 end do
770 end do
771 end if
772 end do
773 end subroutine jacobi_update_lx9
774
775 subroutine jacobi_update_lx8(d, dxt, dyt, dzt, G11, G22, G33, &
776 G12, G13, G23, dfrmd_el, n)
777 integer, parameter :: lx = 8
778 integer, intent(in) :: n
779 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
780 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
781 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
782 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
783 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
784 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
785 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
786 real(kind=rp), intent(in) :: dxt(lx, lx)
787 real(kind=rp), intent(in) :: dyt(lx, lx)
788 real(kind=rp), intent(in) :: dzt(lx, lx)
789 logical, intent(in) :: dfrmd_el(n)
790 integer :: i, j, k, l, e
791
792 d = 0d0
793
794 do e = 1,n
795 do l = 1,lx
796 do k = 1,lx
797 do j = 1,lx
798 do i = 1,lx
799 d(i,j,k,e) = d(i,j,k,e) + &
800 g11(l,j,k,e) * dxt(i,l)**2
801 end do
802 end do
803 end do
804 end do
805 do l = 1,lx
806 do k = 1,lx
807 do j = 1,lx
808 do i = 1,lx
809 d(i,j,k,e) = d(i,j,k,e) + &
810 g22(i,l,k,e) * dyt(j,l)**2
811 end do
812 end do
813 end do
814 end do
815 do l = 1,lx
816 do k = 1,lx
817 do j = 1,lx
818 do i = 1,lx
819 d(i,j,k,e) = d(i,j,k,e) + &
820 g33(i,j,l,e) * dzt(k,l)**2
821 end do
822 end do
823 end do
824 end do
825
826 if (dfrmd_el(e)) then
827 do j = 1,lx,lx-1
828 do k = 1,lx,lx-1
829 d(1,j,k,e) = d(1,j,k,e) &
830 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
831 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
832 d(lx,j,k,e) = d(lx,j,k,e) &
833 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
834 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
835 end do
836 end do
837
838 do i = 1,lx,lx-1
839 do k = 1,lx,lx-1
840 d(i,1,k,e) = d(i,1,k,e) &
841 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
842 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
843 d(i,lx,k,e) = d(i,lx,k,e) &
844 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
845 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
846 end do
847 end do
848 do i = 1,lx,lx-1
849 do j = 1,lx,lx-1
850 d(i,j,1,e) = d(i,j,1,e) &
851 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
852 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
853 d(i,j,lx,e) = d(i,j,lx,e) &
854 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
855 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
856 end do
857 end do
858 end if
859 end do
860 end subroutine jacobi_update_lx8
861
862 subroutine jacobi_update_lx7(d, dxt, dyt, dzt, G11, G22, G33, &
863 G12, G13, G23, dfrmd_el, n)
864 integer, parameter :: lx = 7
865 integer, intent(in) :: n
866 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
867 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
868 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
869 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
870 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
871 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
872 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
873 real(kind=rp), intent(in) :: dxt(lx, lx)
874 real(kind=rp), intent(in) :: dyt(lx, lx)
875 real(kind=rp), intent(in) :: dzt(lx, lx)
876 logical, intent(in) :: dfrmd_el(n)
877 integer :: i, j, k, l, e
878
879 d = 0d0
880
881 do e = 1,n
882 do l = 1,lx
883 do k = 1,lx
884 do j = 1,lx
885 do i = 1,lx
886 d(i,j,k,e) = d(i,j,k,e) + &
887 g11(l,j,k,e) * dxt(i,l)**2
888 end do
889 end do
890 end do
891 end do
892 do l = 1,lx
893 do k = 1,lx
894 do j = 1,lx
895 do i = 1,lx
896 d(i,j,k,e) = d(i,j,k,e) + &
897 g22(i,l,k,e) * dyt(j,l)**2
898 end do
899 end do
900 end do
901 end do
902 do l = 1,lx
903 do k = 1,lx
904 do j = 1,lx
905 do i = 1,lx
906 d(i,j,k,e) = d(i,j,k,e) + &
907 g33(i,j,l,e) * dzt(k,l)**2
908 end do
909 end do
910 end do
911 end do
912
913 if (dfrmd_el(e)) then
914 do j = 1,lx,lx-1
915 do k = 1,lx,lx-1
916 d(1,j,k,e) = d(1,j,k,e) &
917 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
918 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
919 d(lx,j,k,e) = d(lx,j,k,e) &
920 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
921 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
922 end do
923 end do
924
925 do i = 1,lx,lx-1
926 do k = 1,lx,lx-1
927 d(i,1,k,e) = d(i,1,k,e) &
928 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
929 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
930 d(i,lx,k,e) = d(i,lx,k,e) &
931 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
932 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
933 end do
934 end do
935 do i = 1,lx,lx-1
936 do j = 1,lx,lx-1
937 d(i,j,1,e) = d(i,j,1,e) &
938 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
939 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
940 d(i,j,lx,e) = d(i,j,lx,e) &
941 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
942 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
943 end do
944 end do
945 end if
946 end do
947 end subroutine jacobi_update_lx7
948
949 subroutine jacobi_update_lx6(d, dxt, dyt, dzt, G11, G22, G33, &
950 G12, G13, G23, dfrmd_el, n)
951 integer, parameter :: lx = 6
952 integer, intent(in) :: n
953 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
954 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
955 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
956 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
957 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
958 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
959 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
960 real(kind=rp), intent(in) :: dxt(lx, lx)
961 real(kind=rp), intent(in) :: dyt(lx, lx)
962 real(kind=rp), intent(in) :: dzt(lx, lx)
963 logical, intent(in) :: dfrmd_el(n)
964 integer :: i, j, k, l, e
965
966 d = 0d0
967
968 do e = 1,n
969 do l = 1,lx
970 do k = 1,lx
971 do j = 1,lx
972 do i = 1,lx
973 d(i,j,k,e) = d(i,j,k,e) + &
974 g11(l,j,k,e) * dxt(i,l)**2
975 end do
976 end do
977 end do
978 end do
979 do l = 1,lx
980 do k = 1,lx
981 do j = 1,lx
982 do i = 1,lx
983 d(i,j,k,e) = d(i,j,k,e) + &
984 g22(i,l,k,e) * dyt(j,l)**2
985 end do
986 end do
987 end do
988 end do
989 do l = 1,lx
990 do k = 1,lx
991 do j = 1,lx
992 do i = 1,lx
993 d(i,j,k,e) = d(i,j,k,e) + &
994 g33(i,j,l,e) * dzt(k,l)**2
995 end do
996 end do
997 end do
998 end do
999
1000 if (dfrmd_el(e)) then
1001 do j = 1,lx,lx-1
1002 do k = 1,lx,lx-1
1003 d(1,j,k,e) = d(1,j,k,e) &
1004 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1005 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1006 d(lx,j,k,e) = d(lx,j,k,e) &
1007 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1008 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1009 end do
1010 end do
1011
1012 do i = 1,lx,lx-1
1013 do k = 1,lx,lx-1
1014 d(i,1,k,e) = d(i,1,k,e) &
1015 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1016 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1017 d(i,lx,k,e) = d(i,lx,k,e) &
1018 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1019 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1020 end do
1021 end do
1022 do i = 1,lx,lx-1
1023 do j = 1,lx,lx-1
1024 d(i,j,1,e) = d(i,j,1,e) &
1025 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1026 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1027 d(i,j,lx,e) = d(i,j,lx,e) &
1028 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1029 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1030 end do
1031 end do
1032 end if
1033 end do
1034 end subroutine jacobi_update_lx6
1035
1036 subroutine jacobi_update_lx5(d, dxt, dyt, dzt, G11, G22, G33, &
1037 G12, G13, G23, dfrmd_el, n)
1038 integer, parameter :: lx = 5
1039 integer, intent(in) :: n
1040 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1041 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1042 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1043 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1044 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1045 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1046 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1047 real(kind=rp), intent(in) :: dxt(lx, lx)
1048 real(kind=rp), intent(in) :: dyt(lx, lx)
1049 real(kind=rp), intent(in) :: dzt(lx, lx)
1050 logical, intent(in) :: dfrmd_el(n)
1051 integer :: i, j, k, l, e
1052
1053 d = 0d0
1054
1055 do e = 1,n
1056 do l = 1,lx
1057 do k = 1,lx
1058 do j = 1,lx
1059 do i = 1,lx
1060 d(i,j,k,e) = d(i,j,k,e) + &
1061 g11(l,j,k,e) * dxt(i,l)**2
1062 end do
1063 end do
1064 end do
1065 end do
1066 do l = 1,lx
1067 do k = 1,lx
1068 do j = 1,lx
1069 do i = 1,lx
1070 d(i,j,k,e) = d(i,j,k,e) + &
1071 g22(i,l,k,e) * dyt(j,l)**2
1072 end do
1073 end do
1074 end do
1075 end do
1076 do l = 1,lx
1077 do k = 1,lx
1078 do j = 1,lx
1079 do i = 1,lx
1080 d(i,j,k,e) = d(i,j,k,e) + &
1081 g33(i,j,l,e) * dzt(k,l)**2
1082 end do
1083 end do
1084 end do
1085 end do
1086
1087 if (dfrmd_el(e)) then
1088 do j = 1,lx,lx-1
1089 do k = 1,lx,lx-1
1090 d(1,j,k,e) = d(1,j,k,e) &
1091 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1092 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1093 d(lx,j,k,e) = d(lx,j,k,e) &
1094 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1095 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1096 end do
1097 end do
1098
1099 do i = 1,lx,lx-1
1100 do k = 1,lx,lx-1
1101 d(i,1,k,e) = d(i,1,k,e) &
1102 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1103 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1104 d(i,lx,k,e) = d(i,lx,k,e) &
1105 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1106 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1107 end do
1108 end do
1109 do i = 1,lx,lx-1
1110 do j = 1,lx,lx-1
1111 d(i,j,1,e) = d(i,j,1,e) &
1112 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1113 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1114 d(i,j,lx,e) = d(i,j,lx,e) &
1115 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1116 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1117 end do
1118 end do
1119 end if
1120 end do
1121 end subroutine jacobi_update_lx5
1122
1123 subroutine jacobi_update_lx4(d, dxt, dyt, dzt, G11, G22, G33, &
1124 G12, G13, G23, dfrmd_el, n)
1125 integer, parameter :: lx = 4
1126 integer, intent(in) :: n
1127 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1128 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1129 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1130 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1131 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1132 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1133 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1134 real(kind=rp), intent(in) :: dxt(lx, lx)
1135 real(kind=rp), intent(in) :: dyt(lx, lx)
1136 real(kind=rp), intent(in) :: dzt(lx, lx)
1137 logical, intent(in) :: dfrmd_el(n)
1138 integer :: i, j, k, l, e
1139
1140 d = 0d0
1141
1142 do e = 1,n
1143 do l = 1,lx
1144 do k = 1,lx
1145 do j = 1,lx
1146 do i = 1,lx
1147 d(i,j,k,e) = d(i,j,k,e) + &
1148 g11(l,j,k,e) * dxt(i,l)**2
1149 end do
1150 end do
1151 end do
1152 end do
1153 do l = 1,lx
1154 do k = 1,lx
1155 do j = 1,lx
1156 do i = 1,lx
1157 d(i,j,k,e) = d(i,j,k,e) + &
1158 g22(i,l,k,e) * dyt(j,l)**2
1159 end do
1160 end do
1161 end do
1162 end do
1163 do l = 1,lx
1164 do k = 1,lx
1165 do j = 1,lx
1166 do i = 1,lx
1167 d(i,j,k,e) = d(i,j,k,e) + &
1168 g33(i,j,l,e) * dzt(k,l)**2
1169 end do
1170 end do
1171 end do
1172 end do
1173
1174 if (dfrmd_el(e)) then
1175 do j = 1,lx,lx-1
1176 do k = 1,lx,lx-1
1177 d(1,j,k,e) = d(1,j,k,e) &
1178 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1179 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1180 d(lx,j,k,e) = d(lx,j,k,e) &
1181 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1182 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1183 end do
1184 end do
1185
1186 do i = 1,lx,lx-1
1187 do k = 1,lx,lx-1
1188 d(i,1,k,e) = d(i,1,k,e) &
1189 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1190 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1191 d(i,lx,k,e) = d(i,lx,k,e) &
1192 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1193 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1194 end do
1195 end do
1196 do i = 1,lx,lx-1
1197 do j = 1,lx,lx-1
1198 d(i,j,1,e) = d(i,j,1,e) &
1199 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1200 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1201 d(i,j,lx,e) = d(i,j,lx,e) &
1202 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1203 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1204 end do
1205 end do
1206 end if
1207 end do
1208 end subroutine jacobi_update_lx4
1209
1210 subroutine jacobi_update_lx3(d, dxt, dyt, dzt, G11, G22, G33, &
1211 G12, G13, G23, dfrmd_el, n)
1212 integer, parameter :: lx = 3
1213 integer, intent(in) :: n
1214 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1215 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1216 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1217 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1218 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1219 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1220 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1221 real(kind=rp), intent(in) :: dxt(lx, lx)
1222 real(kind=rp), intent(in) :: dyt(lx, lx)
1223 real(kind=rp), intent(in) :: dzt(lx, lx)
1224 logical, intent(in) :: dfrmd_el(n)
1225 integer :: i, j, k, l, e
1226
1227 d = 0d0
1228
1229 do e = 1,n
1230 do l = 1,lx
1231 do k = 1,lx
1232 do j = 1,lx
1233 do i = 1,lx
1234 d(i,j,k,e) = d(i,j,k,e) + &
1235 g11(l,j,k,e) * dxt(i,l)**2
1236 end do
1237 end do
1238 end do
1239 end do
1240 do l = 1,lx
1241 do k = 1,lx
1242 do j = 1,lx
1243 do i = 1,lx
1244 d(i,j,k,e) = d(i,j,k,e) + &
1245 g22(i,l,k,e) * dyt(j,l)**2
1246 end do
1247 end do
1248 end do
1249 end do
1250 do l = 1,lx
1251 do k = 1,lx
1252 do j = 1,lx
1253 do i = 1,lx
1254 d(i,j,k,e) = d(i,j,k,e) + &
1255 g33(i,j,l,e) * dzt(k,l)**2
1256 end do
1257 end do
1258 end do
1259 end do
1260
1261 if (dfrmd_el(e)) then
1262 do j = 1,lx,lx-1
1263 do k = 1,lx,lx-1
1264 d(1,j,k,e) = d(1,j,k,e) &
1265 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1266 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1267 d(lx,j,k,e) = d(lx,j,k,e) &
1268 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1269 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1270 end do
1271 end do
1272
1273 do i = 1,lx,lx-1
1274 do k = 1,lx,lx-1
1275 d(i,1,k,e) = d(i,1,k,e) &
1276 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1277 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1278 d(i,lx,k,e) = d(i,lx,k,e) &
1279 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1280 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1281 end do
1282 end do
1283 do i = 1,lx,lx-1
1284 do j = 1,lx,lx-1
1285 d(i,j,1,e) = d(i,j,1,e) &
1286 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1287 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1288 d(i,j,lx,e) = d(i,j,lx,e) &
1289 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1290 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1291 end do
1292 end do
1293 end if
1294 end do
1295 end subroutine jacobi_update_lx3
1296
1297 subroutine jacobi_update_lx2(d, dxt, dyt, dzt, G11, G22, G33, &
1298 G12, G13, G23, dfrmd_el, n)
1299 integer, parameter :: lx = 2
1300 integer, intent(in) :: n
1301 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1302 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1303 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1304 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1305 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1306 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1307 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1308 real(kind=rp), intent(in) :: dxt(lx, lx)
1309 real(kind=rp), intent(in) :: dyt(lx, lx)
1310 real(kind=rp), intent(in) :: dzt(lx, lx)
1311 logical, intent(in) :: dfrmd_el(n)
1312 integer :: i, j, k, l, e
1313
1314 d = 0d0
1315
1316 do e = 1,n
1317 do l = 1,lx
1318 do k = 1,lx
1319 do j = 1,lx
1320 do i = 1,lx
1321 d(i,j,k,e) = d(i,j,k,e) + &
1322 g11(l,j,k,e) * dxt(i,l)**2
1323 end do
1324 end do
1325 end do
1326 end do
1327 do l = 1,lx
1328 do k = 1,lx
1329 do j = 1,lx
1330 do i = 1,lx
1331 d(i,j,k,e) = d(i,j,k,e) + &
1332 g22(i,l,k,e) * dyt(j,l)**2
1333 end do
1334 end do
1335 end do
1336 end do
1337 do l = 1,lx
1338 do k = 1,lx
1339 do j = 1,lx
1340 do i = 1,lx
1341 d(i,j,k,e) = d(i,j,k,e) + &
1342 g33(i,j,l,e) * dzt(k,l)**2
1343 end do
1344 end do
1345 end do
1346 end do
1347
1348 if (dfrmd_el(e)) then
1349 do j = 1,lx,lx-1
1350 do k = 1,lx,lx-1
1351 d(1,j,k,e) = d(1,j,k,e) &
1352 + g12(1,j,k,e) * dxt(1,1)*dyt(j,j) &
1353 + g13(1,j,k,e) * dxt(1,1)*dzt(k,k)
1354 d(lx,j,k,e) = d(lx,j,k,e) &
1355 + g12(lx,j,k,e) * dxt(lx,lx)*dyt(j,j) &
1356 + g13(lx,j,k,e) * dxt(lx,lx)*dzt(k,k)
1357 end do
1358 end do
1359
1360 do i = 1,lx,lx-1
1361 do k = 1,lx,lx-1
1362 d(i,1,k,e) = d(i,1,k,e) &
1363 + g12(i,1,k,e) * dyt(1,1)*dxt(i,i) &
1364 + g23(i,1,k,e) * dyt(1,1)*dzt(k,k)
1365 d(i,lx,k,e) = d(i,lx,k,e) &
1366 + g12(i,lx,k,e) * dyt(lx,lx)*dxt(i,i) &
1367 + g23(i,lx,k,e) * dyt(lx,lx)*dzt(k,k)
1368 end do
1369 end do
1370 do i = 1,lx,lx-1
1371 do j = 1,lx,lx-1
1372 d(i,j,1,e) = d(i,j,1,e) &
1373 + g13(i,j,1,e) * dzt(1,1)*dxt(i,i) &
1374 + g23(i,j,1,e) * dzt(1,1)*dyt(j,j)
1375 d(i,j,lx,e) = d(i,j,lx,e) &
1376 + g13(i,j,lx,e) * dzt(lx,lx)*dxt(i,i) &
1377 + g23(i,j,lx,e) * dzt(lx,lx)*dyt(j,j)
1378 end do
1379 end do
1380 end if
1381 end do
1382 end subroutine jacobi_update_lx2
1383
1384end 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:124
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:777
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:864
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:516
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:603
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:169
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:342
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:429
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:1038
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:690
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:255
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:1125
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:1212
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:951
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:1299
math
Definition math.f90:60
math::addcol3
subroutine, public addcol3(a, b, c, n)
Returns .
Definition math.f90:958
math::invcol1
subroutine, public invcol1(a, n)
Invert a vector .
Definition math.f90:614
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:854
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:867
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:55
dofmap::dofmap_t
Definition dofmap.f90:53
gather_scatter::gs_t
Definition gather_scatter.f90:68
jacobi::jacobi_t
Defines a jacobi preconditioner.
Definition pc_jacobi.f90:45
precon::pc_t
Defines a canonical Krylov preconditioner.
Definition precon.f90:40