Neko  0.9.99
A portable framework for high-order spectral element flow simulations
pipecg.f90
Go to the documentation of this file.
1 ! Copyright (c) 2021-2024, 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 pipecg
36  use precon, only : pc_t
37  use ax_product, only : ax_t
38  use num_types, only: rp
39  use field, only : field_t
40  use coefs, only : coef_t
41  use gather_scatter, only : gs_t, gs_op_add
42  use bc, only : bc_list_t, bc_list_apply
43  use math, only : glsc3, rzero, copy, abscmp
44  use comm
45  implicit none
46  private
47 
48  integer, parameter :: pipecg_p_space = 7
49 
51  type, public, extends(ksp_t) :: pipecg_t
52  real(kind=rp), allocatable :: p(:)
53  real(kind=rp), allocatable :: q(:)
54  real(kind=rp), allocatable :: r(:)
55  real(kind=rp), allocatable :: s(:)
56  real(kind=rp), allocatable :: u(:,:)
57  real(kind=rp), allocatable :: w(:)
58  real(kind=rp), allocatable :: z(:)
59  real(kind=rp), allocatable :: mi(:)
60  real(kind=rp), allocatable :: ni(:)
61  contains
63  procedure, pass(this) :: init => pipecg_init
65  procedure, pass(this) :: free => pipecg_free
67  procedure, pass(this) :: solve => pipecg_solve
69  procedure, pass(this) :: solve_coupled => pipecg_solve_coupled
70  end type pipecg_t
71 
72 contains
73 
75  subroutine pipecg_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
76  class(pipecg_t), intent(inout) :: this
77  integer, intent(in) :: max_iter
78  class(pc_t), optional, intent(in), target :: M
79  integer, intent(in) :: n
80  real(kind=rp), optional, intent(in) :: rel_tol
81  real(kind=rp), optional, intent(in) :: abs_tol
82  logical, optional, intent(in) :: monitor
83 
84  call this%free()
85 
86  allocate(this%p(n))
87  allocate(this%q(n))
88  allocate(this%r(n))
89  allocate(this%s(n))
90  allocate(this%u(n,pipecg_p_space+1))
91  allocate(this%w(n))
92  allocate(this%z(n))
93  allocate(this%mi(n))
94  allocate(this%ni(n))
95  if (present(m)) then
96  this%M => m
97  end if
98 
99  if (present(rel_tol) .and. present(abs_tol) .and. present(monitor)) then
100  call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
101  else if (present(rel_tol) .and. present(abs_tol)) then
102  call this%ksp_init(max_iter, rel_tol, abs_tol)
103  else if (present(monitor) .and. present(abs_tol)) then
104  call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
105  else if (present(rel_tol) .and. present(monitor)) then
106  call this%ksp_init(max_iter, rel_tol, monitor = monitor)
107  else if (present(rel_tol)) then
108  call this%ksp_init(max_iter, rel_tol = rel_tol)
109  else if (present(abs_tol)) then
110  call this%ksp_init(max_iter, abs_tol = abs_tol)
111  else if (present(monitor)) then
112  call this%ksp_init(max_iter, monitor = monitor)
113  else
114  call this%ksp_init(max_iter)
115  end if
116 
117  end subroutine pipecg_init
118 
120  subroutine pipecg_free(this)
121  class(pipecg_t), intent(inout) :: this
122 
123  call this%ksp_free()
124 
125  if (allocated(this%p)) then
126  deallocate(this%p)
127  end if
128  if (allocated(this%q)) then
129  deallocate(this%q)
130  end if
131  if (allocated(this%r)) then
132  deallocate(this%r)
133  end if
134  if (allocated(this%s)) then
135  deallocate(this%s)
136  end if
137  if (allocated(this%u)) then
138  deallocate(this%u)
139  end if
140  if (allocated(this%w)) then
141  deallocate(this%w)
142  end if
143  if (allocated(this%z)) then
144  deallocate(this%z)
145  end if
146  if (allocated(this%mi)) then
147  deallocate(this%mi)
148  end if
149  if (allocated(this%ni)) then
150  deallocate(this%ni)
151  end if
152 
153  nullify(this%M)
154 
155 
156  end subroutine pipecg_free
157 
159  function pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_results)
160  class(pipecg_t), intent(inout) :: this
161  class(ax_t), intent(in) :: ax
162  type(field_t), intent(inout) :: x
163  integer, intent(in) :: n
164  real(kind=rp), dimension(n), intent(in) :: f
165  type(coef_t), intent(inout) :: coef
166  type(bc_list_t), intent(in) :: blst
167  type(gs_t), intent(inout) :: gs_h
168  type(ksp_monitor_t) :: ksp_results
169  integer, optional, intent(in) :: niter
170  integer :: iter, max_iter, i, j, k, ierr, p_cur, p_prev, u_prev
171  real(kind=rp) :: rnorm, rtr, reduction(3), norm_fac
172  real(kind=rp) :: alpha(pipecg_p_space), beta(pipecg_p_space)
173  real(kind=rp) :: gamma1, gamma2, delta
174  real(kind=rp) :: tmp1, tmp2, tmp3, x_plus(neko_blk_size)
175  type(mpi_request) :: request
176  type(mpi_status) :: status
177 
178  if (present(niter)) then
179  max_iter = niter
180  else
181  max_iter = this%max_iter
182  end if
183  norm_fac = 1.0_rp / sqrt(coef%volume)
184 
185  associate(p => this%p, q => this%q, r => this%r, s => this%s, &
186  u => this%u, w => this%w, z => this%z, mi => this%mi, ni => this%ni)
187 
188  p_prev = pipecg_p_space
189  u_prev = pipecg_p_space+1
190  p_cur = 1
191  call rzero(x%x, n)
192  call rzero(z, n)
193  call rzero(q, n)
194  call rzero(p, n)
195  call rzero(s, n)
196  call copy(r, f, n)
197  call this%M%solve(u(1,u_prev), r, n)
198  call ax%compute(w, u(1,u_prev), coef, x%msh, x%Xh)
199  call gs_h%op(w, n, gs_op_add)
200  call bc_list_apply(blst, w, n)
201 
202  rtr = glsc3(r, coef%mult, r, n)
203  rnorm = sqrt(rtr)*norm_fac
204  ksp_results%res_start = rnorm
205  ksp_results%res_final = rnorm
206  ksp_results%iter = 0
207  if(abscmp(rnorm, 0.0_rp)) return
208 
209  gamma1 = 0.0_rp
210  tmp1 = 0.0_rp
211  tmp2 = 0.0_rp
212  tmp3 = 0.0_rp
213  do i = 1, n
214  tmp1 = tmp1 + r(i) * coef%mult(i,1,1,1) * u(i,u_prev)
215  tmp2 = tmp2 + w(i) * coef%mult(i,1,1,1) * u(i,u_prev)
216  tmp3 = tmp3 + r(i) * coef%mult(i,1,1,1) * r(i)
217  end do
218  reduction(1) = tmp1
219  reduction(2) = tmp2
220  reduction(3) = tmp3
221 
222  call this%monitor_start('PipeCG')
223  do iter = 1, max_iter
224  call mpi_iallreduce(mpi_in_place, reduction, 3, &
225  mpi_real_precision, mpi_sum, neko_comm, request, ierr)
226 
227  call this%M%solve(mi, w, n)
228  call ax%compute(ni, mi, coef, x%msh, x%Xh)
229  call gs_h%op(ni, n, gs_op_add)
230  call bc_list_apply(blst, ni, n)
231 
232  call mpi_wait(request, status, ierr)
233  gamma2 = gamma1
234  gamma1 = reduction(1)
235  delta = reduction(2)
236  rtr = reduction(3)
237 
238  rnorm = sqrt(rtr)*norm_fac
239  call this%monitor_iter(iter, rnorm)
240  if (rnorm .lt. this%abs_tol) exit
241 
242  if (iter .gt. 1) then
243  beta(p_cur) = gamma1 / gamma2
244  alpha(p_cur) = gamma1 / (delta - (beta(p_cur) * gamma1/alpha(p_prev)))
245  else
246  beta(p_cur) = 0.0_rp
247  alpha(p_cur) = gamma1/delta
248  end if
249 
250  tmp1 = 0.0_rp
251  tmp2 = 0.0_rp
252  tmp3 = 0.0_rp
253  do i = 0, n, neko_blk_size
254  if (i + neko_blk_size .le. n) then
255  do k = 1, neko_blk_size
256  z(i+k) = beta(p_cur) * z(i+k) + ni(i+k)
257  q(i+k) = beta(p_cur) * q(i+k) + mi(i+k)
258  s(i+k) = beta(p_cur) * s(i+k) + w(i+k)
259  r(i+k) = r(i+k) - alpha(p_cur) * s(i+k)
260  u(i+k,p_cur) = u(i+k,u_prev) - alpha(p_cur) * q(i+k)
261  w(i+k) = w(i+k) - alpha(p_cur) * z(i+k)
262  tmp1 = tmp1 + r(i+k) * coef%mult(i+k,1,1,1) * u(i+k,p_cur)
263  tmp2 = tmp2 + w(i+k) * coef%mult(i+k,1,1,1) * u(i+k,p_cur)
264  tmp3 = tmp3 + r(i+k) * coef%mult(i+k,1,1,1) * r(i+k)
265  end do
266  else
267  do k = 1, n-i
268  z(i+k) = beta(p_cur) * z(i+k) + ni(i+k)
269  q(i+k) = beta(p_cur) * q(i+k) + mi(i+k)
270  s(i+k) = beta(p_cur) * s(i+k) + w(i+k)
271  r(i+k) = r(i+k) - alpha(p_cur) * s(i+k)
272  u(i+k,p_cur) = u(i+k,u_prev) - alpha(p_cur) * q(i+k)
273  w(i+k) = w(i+k) - alpha(p_cur) * z(i+k)
274  tmp1 = tmp1 + r(i+k) * coef%mult(i+k,1,1,1) * u(i+k,p_cur)
275  tmp2 = tmp2 + w(i+k) * coef%mult(i+k,1,1,1) * u(i+k,p_cur)
276  tmp3 = tmp3 + r(i+k) * coef%mult(i+k,1,1,1) * r(i+k)
277  end do
278  end if
279  end do
280 
281  reduction(1) = tmp1
282  reduction(2) = tmp2
283  reduction(3) = tmp3
284 
285  if (p_cur .eq. pipecg_p_space) then
286  do i = 0, n, neko_blk_size
287  if (i + neko_blk_size .le. n) then
288  do k = 1, neko_blk_size
289  x_plus(k) = 0.0_rp
290  end do
291  p_prev = pipecg_p_space+1
292  do j = 1, p_cur
293  do k = 1, neko_blk_size
294  p(i+k) = beta(j) * p(i+k) + u(i+k,p_prev)
295  x_plus(k) = x_plus(k) + alpha(j) * p(i+k)
296  end do
297  p_prev = j
298  end do
299  do k = 1, neko_blk_size
300  x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
301  u(i+k,pipecg_p_space+1) = u(i+k,pipecg_p_space)
302  end do
303  else
304  do k = 1, n-i
305  x_plus(1) = 0.0_rp
306  p_prev = pipecg_p_space + 1
307  do j = 1, p_cur
308  p(i+k) = beta(j) * p(i+k) + u(i+k,p_prev)
309  x_plus(1) = x_plus(1) + alpha(j) * p(i+k)
310  p_prev = j
311  end do
312  x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
313  u(i+k,pipecg_p_space+1) = u(i+k,pipecg_p_space)
314  end do
315  end if
316  end do
317  p_prev = p_cur
318  u_prev = pipecg_p_space+1
319  alpha(1) = alpha(p_cur)
320  beta(1) = beta(p_cur)
321  p_cur = 1
322  else
323  u_prev = p_cur
324  p_prev = p_cur
325  p_cur = p_cur + 1
326  end if
327  end do
328 
329  if ( p_cur .ne. 1) then
330  do i = 0, n, neko_blk_size
331  if (i + neko_blk_size .le. n) then
332  do k = 1, neko_blk_size
333  x_plus(k) = 0.0_rp
334  end do
335  p_prev = pipecg_p_space+1
336  do j = 1, p_cur
337  do k = 1, neko_blk_size
338  p(i+k) = beta(j) * p(i+k) + u(i+k,p_prev)
339  x_plus(k) = x_plus(k) + alpha(j) * p(i+k)
340  end do
341  p_prev = j
342  end do
343  do k = 1, neko_blk_size
344  x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
345  u(i+k,pipecg_p_space+1) = u(i+k,pipecg_p_space)
346  end do
347  else
348  do k = 1, n-i
349  x_plus(1) = 0.0_rp
350  p_prev = pipecg_p_space + 1
351  do j = 1, p_cur
352  p(i+k) = beta(j) * p(i+k) + u(i+k,p_prev)
353  x_plus(1) = x_plus(1) + alpha(j) * p(i+k)
354  p_prev = j
355  end do
356  x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
357  u(i+k,pipecg_p_space+1) = u(i+k,pipecg_p_space)
358  end do
359  end if
360  end do
361  end if
362  call this%monitor_stop()
363  ksp_results%res_final = rnorm
364  ksp_results%iter = iter
365 
366  end associate
367 
368  end function pipecg_solve
369 
371  function pipecg_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
372  n, coef, blstx, blsty, blstz, gs_h, niter) result(ksp_results)
373  class(pipecg_t), intent(inout) :: this
374  class(ax_t), intent(in) :: ax
375  type(field_t), intent(inout) :: x
376  type(field_t), intent(inout) :: y
377  type(field_t), intent(inout) :: z
378  integer, intent(in) :: n
379  real(kind=rp), dimension(n), intent(in) :: fx
380  real(kind=rp), dimension(n), intent(in) :: fy
381  real(kind=rp), dimension(n), intent(in) :: fz
382  type(coef_t), intent(inout) :: coef
383  type(bc_list_t), intent(in) :: blstx
384  type(bc_list_t), intent(in) :: blsty
385  type(bc_list_t), intent(in) :: blstz
386  type(gs_t), intent(inout) :: gs_h
387  type(ksp_monitor_t), dimension(3) :: ksp_results
388  integer, optional, intent(in) :: niter
389 
390  ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
391  ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
392  ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
393 
394  end function pipecg_solve_coupled
395 
396 end module pipecg
397 
398 
Defines a Matrix-vector product.
Definition: ax.f90:34
Defines a boundary condition.
Definition: bc.f90:34
Coefficients.
Definition: coef.f90:34
Definition: comm.F90:1
type(mpi_comm) neko_comm
MPI communicator.
Definition: comm.F90:16
type(mpi_datatype) mpi_real_precision
MPI type for working precision of REAL types.
Definition: comm.F90:23
Defines a field.
Definition: field.f90:34
Gather-scatter.
Implements the base abstract type for Krylov solvers plus helper types.
Definition: krylov.f90:34
integer, parameter, public ksp_max_iter
Maximum number of iters.
Definition: krylov.f90:51
Definition: math.f90:60
real(kind=rp) function, public glsc3(a, b, c, n)
Weighted inner product .
Definition: math.f90:895
subroutine, public copy(a, b, n)
Copy a vector .
Definition: math.f90:239
subroutine, public rzero(a, n)
Zero a real vector.
Definition: math.f90:195
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Defines a pipelined Conjugate Gradient methods.
Definition: pipecg.f90:34
subroutine pipecg_free(this)
Deallocate a pipelined PCG solver.
Definition: pipecg.f90:121
integer, parameter pipecg_p_space
Definition: pipecg.f90:48
subroutine pipecg_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
Initialise a pipelined PCG solver.
Definition: pipecg.f90:76
type(ksp_monitor_t) function, dimension(3) pipecg_solve_coupled(this, Ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
Pipelined PCG coupled solve.
Definition: pipecg.f90:373
type(ksp_monitor_t) function pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Pipelined PCG solve.
Definition: pipecg.f90:160
Krylov preconditioner.
Definition: precon.f90:34
Base type for a matrix-vector product providing .
Definition: ax.f90:43
A list of boundary conditions.
Definition: bc.f90:104
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
Type for storing initial and final residuals in a Krylov solver.
Definition: krylov.f90:56
Base abstract type for a canonical Krylov method, solving .
Definition: krylov.f90:66
Pipelined preconditioned conjugate gradient method.
Definition: pipecg.f90:51
Defines a canonical Krylov preconditioner.
Definition: precon.f90:40