Neko  0.9.0
A portable framework for high-order spectral element flow simulations
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 !
34 module jacobi
35  use math
36  use precon
37  use coefs, only : coef_t
38  use num_types, only : rp
39  use dofmap
40  use gather_scatter
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 
57 contains
58 
59  subroutine jacobi_init(this, coef, dof, gs_h)
60  class(jacobi_t), intent(inout) :: this
61  type(coef_t), intent(inout), target :: coef
62  type(dofmap_t), intent(inout), 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 
1384 end module jacobi
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_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_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_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_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_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_init
subroutine jacobi_init(this, coef, dof, gs_h)
Definition: pc_jacobi.f90:60
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
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_free
subroutine jacobi_free(this)
Definition: pc_jacobi.f90:75
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_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_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_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_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_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_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_lx11
subroutine jacobi_update_lx11(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, dfrmd_el, n)
Definition: pc_jacobi.f90:516
math
Definition: math.f90:60
math::addcol3
subroutine, public addcol3(a, b, c, n)
Returns .
Definition: math.f90:801
math::invcol1
subroutine, public invcol1(a, n)
Invert a vector .
Definition: math.f90:475
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition: math.f90:729
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition: math.f90:742
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:58
jacobi::jacobi_t
Defines a jacobi preconditioner.
Definition: pc_jacobi.f90:45
precon::pc_t
Defines a canonical Krylov preconditioner.
Definition: precon.f90:40