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(:)
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
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)
114 call this%ksp_init(max_iter)
121 class(
pipecg_t),
intent(inout) :: this
125 if (
allocated(this%p))
then
128 if (
allocated(this%q))
then
131 if (
allocated(this%r))
then
134 if (
allocated(this%s))
then
137 if (
allocated(this%u))
then
140 if (
allocated(this%w))
then
143 if (
allocated(this%z))
then
146 if (
allocated(this%mi))
then
149 if (
allocated(this%ni))
then
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
167 type(
gs_t),
intent(inout) :: gs_h
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
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
178 if (
present(niter))
then
181 max_iter = this%max_iter
183 norm_fac = 1.0_rp / sqrt(coef%volume)
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)
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)
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
207 if(
abscmp(rnorm, 0.0_rp))
return
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)
222 call this%monitor_start(
'PipeCG')
223 do iter = 1, max_iter
224 call mpi_iallreduce(mpi_in_place, reduction, 3, &
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)
232 call mpi_wait(request, status, ierr)
234 gamma1 = reduction(1)
238 rnorm = sqrt(rtr)*norm_fac
239 call this%monitor_iter(iter, rnorm)
240 if (rnorm .lt. this%abs_tol)
exit
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)))
247 alpha(p_cur) = gamma1/delta
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)
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)
286 do i = 0, n, neko_blk_size
287 if (i + neko_blk_size .le. n)
then
288 do k = 1, neko_blk_size
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)
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)
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)
312 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
319 alpha(1) = alpha(p_cur)
320 beta(1) = beta(p_cur)
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
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)
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)
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)
356 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
362 call this%monitor_stop()
363 ksp_results%res_final = rnorm
364 ksp_results%iter = iter
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
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)
Defines a Matrix-vector product.
Defines a boundary condition.
type(mpi_comm) neko_comm
MPI communicator.
type(mpi_datatype) mpi_real_precision
MPI type for working precision of REAL types.
Implements the base abstract type for Krylov solvers plus helper types.
integer, parameter, public ksp_max_iter
Maximum number of iters.
real(kind=rp) function, public glsc3(a, b, c, n)
Weighted inner product .
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
integer, parameter, public rp
Global precision used in computations.
Defines a pipelined Conjugate Gradient methods.
subroutine pipecg_free(this)
Deallocate a pipelined PCG solver.
integer, parameter pipecg_p_space
subroutine pipecg_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
Initialise a pipelined PCG solver.
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.
type(ksp_monitor_t) function pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Pipelined PCG solve.
Base type for a matrix-vector product providing .
A list of boundary conditions.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Type for storing initial and final residuals in a Krylov solver.
Base abstract type for a canonical Krylov method, solving .
Pipelined preconditioned conjugate gradient method.
Defines a canonical Krylov preconditioner.