52 real(kind=
rp),
allocatable :: w(:)
53 real(kind=
rp),
allocatable :: r(:)
54 real(kind=
rp),
allocatable :: p(:,:)
55 real(kind=
rp),
allocatable :: z(:)
56 real(kind=
rp),
allocatable :: alpha(:)
67 subroutine cg_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
68 class(
cg_t),
intent(inout),
target :: this
69 integer,
intent(in) :: max_iter
70 class(
pc_t),
optional,
intent(in),
target :: M
71 integer,
intent(in) :: n
72 real(kind=
rp),
optional,
intent(in) :: rel_tol
73 real(kind=
rp),
optional,
intent(in) :: abs_tol
74 logical,
optional,
intent(in) :: monitor
87 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
88 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
89 else if (
present(rel_tol) .and.
present(abs_tol))
then
90 call this%ksp_init(max_iter, rel_tol, abs_tol)
91 else if (
present(monitor) .and.
present(abs_tol))
then
92 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
93 else if (
present(rel_tol) .and.
present(monitor))
then
94 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
95 else if (
present(rel_tol))
then
96 call this%ksp_init(max_iter, rel_tol = rel_tol)
97 else if (
present(abs_tol))
then
98 call this%ksp_init(max_iter, abs_tol = abs_tol)
99 else if (
present(monitor))
then
100 call this%ksp_init(max_iter, monitor = monitor)
102 call this%ksp_init(max_iter)
109 class(
cg_t),
intent(inout) :: this
113 if (
allocated(this%w))
then
117 if (
allocated(this%r))
then
121 if (
allocated(this%p))
then
125 if (
allocated(this%z))
then
129 if (
allocated(this%alpha))
then
130 deallocate(this%alpha)
138 function cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
139 class(
cg_t),
intent(inout) :: this
140 class(
ax_t),
intent(in) :: ax
141 type(
field_t),
intent(inout) :: x
142 integer,
intent(in) :: n
143 real(kind=
rp),
dimension(n),
intent(in) :: f
144 type(
coef_t),
intent(inout) :: coef
146 type(
gs_t),
intent(inout) :: gs_h
148 integer,
optional,
intent(in) :: niter
149 integer :: iter, max_iter, i, j, k, p_cur, p_prev
150 real(kind=
rp) :: rnorm, rtr, rtz2, rtz1, x_plus(neko_blk_size)
151 real(kind=
rp) :: beta, pap, norm_fac
153 if (
present(niter))
then
156 max_iter = this%max_iter
158 norm_fac = 1.0_rp / sqrt(coef%volume)
160 associate(w => this%w, r => this%r, p => this%p, &
161 z => this%z, alpha => this%alpha)
168 rtr =
glsc3(r, coef%mult, r, n)
169 rnorm = sqrt(rtr) * norm_fac
170 ksp_results%res_start = rnorm
171 ksp_results%res_final = rnorm
175 if(
abscmp(rnorm, 0.0_rp))
return
176 call this%monitor_start(
'CG')
177 do iter = 1, max_iter
178 call this%M%solve(z, r, n)
180 rtz1 =
glsc3(r, coef%mult, z, n)
183 if (iter .eq. 1) beta = 0.0_rp
185 p(i,p_cur) = z(i) + beta * p(i,p_prev)
188 call ax%compute(w, p(1,p_cur), coef, x%msh, x%Xh)
189 call gs_h%op(w, n, gs_op_add)
192 pap =
glsc3(w, coef%mult, p(1,p_cur), n)
194 alpha(p_cur) = rtz1 / pap
196 rnorm = sqrt(rtr) * norm_fac
197 call this%monitor_iter(iter, rnorm)
200 (rnorm .lt. this%abs_tol) .or. iter .eq. max_iter)
then
201 do i = 0, n, neko_blk_size
202 if (i + neko_blk_size .le. n)
then
203 do k = 1, neko_blk_size
207 do k = 1, neko_blk_size
208 x_plus(k) = x_plus(k) + alpha(j) * p(i+k,j)
211 do k = 1, neko_blk_size
212 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
218 x_plus(1) = x_plus(1) + alpha(j) * p(i+k,j)
220 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
226 if (rnorm .lt. this%abs_tol)
exit
233 call this%monitor_stop()
234 ksp_results%res_final = rnorm
235 ksp_results%iter = iter
240 integer,
intent(in) :: n
241 real(kind=rp),
intent(inout) :: r(n), rtr
243 real(kind=rp),
intent(in) ::mult(n), w(n), alpha
248 r(i) = r(i) - alpha*w(i)
249 tmp = tmp + r(i) * r(i) * mult(i)
251 call mpi_allreduce(mpi_in_place, tmp, 1, &
252 mpi_extra_precision, mpi_sum, neko_comm, ierr)
259 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
260 class(
cg_t),
intent(inout) :: this
261 class(ax_t),
intent(in) :: ax
262 type(field_t),
intent(inout) :: x
263 type(field_t),
intent(inout) :: y
264 type(field_t),
intent(inout) :: z
265 integer,
intent(in) :: n
266 real(kind=rp),
dimension(n),
intent(in) :: fx
267 real(kind=rp),
dimension(n),
intent(in) :: fy
268 real(kind=rp),
dimension(n),
intent(in) :: fz
269 type(coef_t),
intent(inout) :: coef
270 type(bc_list_t),
intent(in) :: blstx
271 type(bc_list_t),
intent(in) :: blsty
272 type(bc_list_t),
intent(in) :: blstz
273 type(gs_t),
intent(inout) :: gs_h
274 type(ksp_monitor_t),
dimension(3) :: ksp_results
275 integer,
optional,
intent(in) :: niter
277 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
278 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
279 ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
Defines a Matrix-vector product.
Defines a boundary condition.
Defines various Conjugate Gradient methods.
subroutine cg_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
Initialise a standard PCG solver.
subroutine cg_free(this)
Deallocate a standard PCG solver.
type(ksp_monitor_t) function, dimension(3) cg_solve_coupled(this, Ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
Standard PCG coupled solve.
integer, parameter cg_p_space
type(ksp_monitor_t) function cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Standard PCG solve.
subroutine second_cg_part(rtr, r, mult, w, alpha, n)
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 xp
integer, parameter, public rp
Global precision used in computations.
Base type for a matrix-vector product providing .
A list of boundary conditions.
Standard 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.