50 real(kind=
rp),
allocatable :: w1(:)
51 real(kind=
rp),
allocatable :: w2(:)
52 real(kind=
rp),
allocatable :: w3(:)
53 real(kind=
rp),
allocatable :: r1(:)
54 real(kind=
rp),
allocatable :: r2(:)
55 real(kind=
rp),
allocatable :: r3(:)
56 real(kind=
rp),
allocatable :: p1(:)
57 real(kind=
rp),
allocatable :: p2(:)
58 real(kind=
rp),
allocatable :: p3(:)
59 real(kind=
rp),
allocatable :: z1(:)
60 real(kind=
rp),
allocatable :: z2(:)
61 real(kind=
rp),
allocatable :: z3(:)
62 real(kind=
rp),
allocatable :: tmp(:)
73 subroutine cg_cpld_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
75 integer,
intent(in) :: max_iter
76 class(
pc_t),
optional,
intent(in),
target :: M
77 integer,
intent(in) :: n
78 real(kind=
rp),
optional,
intent(in) :: rel_tol
79 real(kind=
rp),
optional,
intent(in) :: abs_tol
80 logical,
optional,
intent(in) :: monitor
102 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
103 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
104 else if (
present(rel_tol) .and.
present(abs_tol))
then
105 call this%ksp_init(max_iter, rel_tol, abs_tol)
106 else if (
present(monitor) .and.
present(abs_tol))
then
107 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
108 else if (
present(rel_tol) .and.
present(monitor))
then
109 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
110 else if (
present(rel_tol))
then
111 call this%ksp_init(max_iter, rel_tol = rel_tol)
112 else if (
present(abs_tol))
then
113 call this%ksp_init(max_iter, abs_tol = abs_tol)
114 else if (
present(monitor))
then
115 call this%ksp_init(max_iter, monitor = monitor)
117 call this%ksp_init(max_iter)
128 if (
allocated(this%w1))
then
132 if (
allocated(this%w2))
then
136 if (
allocated(this%w3))
then
140 if (
allocated(this%r1))
then
144 if (
allocated(this%r2))
then
148 if (
allocated(this%r3))
then
152 if (
allocated(this%p1))
then
156 if (
allocated(this%p2))
then
160 if (
allocated(this%p3))
then
164 if (
allocated(this%z1))
then
168 if (
allocated(this%z2))
then
172 if (
allocated(this%z3))
then
176 if (
allocated(this%tmp))
then
184 function cg_cpld_nop(this, Ax, x, f, n, coef, blst, gs_h, niter) &
187 class(
ax_t),
intent(in) :: ax
188 type(
field_t),
intent(inout) :: x
189 integer,
intent(in) :: n
190 real(kind=
rp),
dimension(n),
intent(in) :: f
191 type(
coef_t),
intent(inout) :: coef
193 type(
gs_t),
intent(inout) :: gs_h
195 integer,
optional,
intent(in) :: niter
198 call neko_error(
'The cpldcg solver is only defined for coupled solves')
200 ksp_results%res_final = 0.0
206 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
208 class(
ax_t),
intent(in) :: ax
209 type(
field_t),
intent(inout) :: x
210 type(
field_t),
intent(inout) :: y
211 type(
field_t),
intent(inout) :: z
212 integer,
intent(in) :: n
213 real(kind=
rp),
dimension(n),
intent(in) :: fx
214 real(kind=
rp),
dimension(n),
intent(in) :: fy
215 real(kind=
rp),
dimension(n),
intent(in) :: fz
216 type(
coef_t),
intent(inout) :: coef
220 type(
gs_t),
intent(inout) :: gs_h
222 integer,
optional,
intent(in) :: niter
223 real(kind=
rp),
parameter :: one = 1.0
224 real(kind=
rp),
parameter :: zero = 0.0
225 integer :: i, iter, max_iter
226 real(kind=
rp) :: rnorm, rtr, rtr0, rtz2, rtz1
227 real(kind=
rp) :: beta, pap, alpha, alphm, norm_fac
229 if (
present(niter))
then
232 max_iter = this%max_iter
234 norm_fac = one / coef%volume
236 associate(p1 => this%p1, p2 => this%p2, p3 => this%p3, z1 => this%z1, &
237 z2 => this%z2, z3 => this%z3, r1 => this%r1, r2 => this%r2, &
238 r3 => this%r3, tmp => this%tmp, w1 => this%w1, w2 => this%w2, &
242 do concurrent(i = 1:n)
243 x%x(i,1,1,1) = 0.0_rp
244 y%x(i,1,1,1) = 0.0_rp
245 z%x(i,1,1,1) = 0.0_rp
255 tmp(i) = r1(i)**2 + r2(i)**2 + r3(i)**2
258 rtr =
glsc3(tmp, coef%mult, coef%binv, n)
259 rnorm = sqrt(rtr*norm_fac)
260 ksp_results%res_start = rnorm
261 ksp_results%res_final = rnorm
263 if (rnorm .eq. zero)
return
265 call this%monitor_start(
'cpldCG')
266 do iter = 1, max_iter
267 call this%M%solve(z1, this%r1, n)
268 call this%M%solve(z2, this%r2, n)
269 call this%M%solve(z3, this%r3, n)
272 do concurrent(i = 1:n)
273 this%tmp(i) = z1(i) * r1(i) &
278 rtz1 =
glsc2(tmp, coef%mult, n)
281 if (iter .eq. 1) beta = zero
282 do concurrent(i = 1:n)
283 p1(i) = p1(i) * beta + z1(i)
284 p2(i) = p2(i) * beta + z2(i)
285 p3(i) = p3(i) * beta + z3(i)
288 call ax%compute_vector(w1, w2, w3, p1, p2, p3, coef, x%msh, x%Xh)
289 call gs_h%op(w1, n, gs_op_add)
290 call gs_h%op(w2, n, gs_op_add)
291 call gs_h%op(w3, n, gs_op_add)
296 do concurrent(i = 1:n)
297 tmp(i) = w1(i) * p1(i) &
302 pap =
glsc2(tmp, coef%mult, n)
306 do concurrent(i = 1:n)
307 x%x(i,1,1,1) = x%x(i,1,1,1) + alpha * p1(i)
308 y%x(i,1,1,1) = y%x(i,1,1,1) + alpha * p2(i)
309 z%x(i,1,1,1) = z%x(i,1,1,1) + alpha * p3(i)
310 r1(i) = r1(i) + alphm * w1(i)
311 r2(i) = r2(i) + alphm * w2(i)
312 r3(i) = r3(i) + alphm * w3(i)
313 tmp(i) = r1(i)**2 + r2(i)**2 + r3(i)**2
316 rtr =
glsc3(tmp, coef%mult, coef%binv, n)
317 if (iter .eq. 1) rtr0 = rtr
318 rnorm = sqrt(rtr * norm_fac)
319 call this%monitor_iter(iter, rnorm)
320 if (rnorm .lt. this%abs_tol)
then
325 call this%monitor_stop()
326 ksp_results%res_final = rnorm
327 ksp_results%iter = iter
Defines a Matrix-vector product.
Defines a boundary condition.
Defines a coupled Conjugate Gradient methods.
type(ksp_monitor_t) function, dimension(3) cg_cpld_solve(this, Ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
Coupled PCG solve.
subroutine cg_cpld_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
Initialise a coupled PCG solver.
type(ksp_monitor_t) function cg_cpld_nop(this, Ax, x, f, n, coef, blst, gs_h, niter)
subroutine cg_cpld_free(this)
Deallocate a coupled PCG solver.
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 .
real(kind=rp) function, public glsc2(a, b, n)
Weighted inner product .
integer, parameter, public rp
Global precision used in computations.
Base type for a matrix-vector product providing .
A list of boundary conditions.
Coupled preconditioned conjugate gradient method.
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 .
Defines a canonical Krylov preconditioner.