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