52 real(kind=
rp),
allocatable :: r(:)
53 real(kind=
rp),
allocatable :: p(:)
54 real(kind=
rp),
allocatable :: pr(:,:)
66 subroutine cacg_init(this, n, max_iter, M, s, rel_tol, abs_tol, monitor)
67 class(
cacg_t),
intent(inout) :: this
68 class(
pc_t),
optional,
intent(inout),
target :: M
69 integer,
intent(in) :: n
70 integer,
intent(in) :: max_iter
71 real(kind=
rp),
optional,
intent(inout) :: rel_tol
72 real(kind=
rp),
optional,
intent(inout) :: abs_tol
73 logical,
optional,
intent(in) :: monitor
74 integer,
optional,
intent(inout) :: s
84 & be aware of potential instabilities")
89 allocate(this%PR(n,4*this%s+1))
94 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
95 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
96 else if (
present(rel_tol) .and.
present(abs_tol))
then
97 call this%ksp_init(max_iter, rel_tol, abs_tol)
98 else if (
present(monitor) .and.
present(abs_tol))
then
99 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
100 else if (
present(rel_tol) .and.
present(monitor))
then
101 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
102 else if (
present(rel_tol))
then
103 call this%ksp_init(max_iter, rel_tol = rel_tol)
104 else if (
present(abs_tol))
then
105 call this%ksp_init(max_iter, abs_tol = abs_tol)
106 else if (
present(monitor))
then
107 call this%ksp_init(max_iter, monitor = monitor)
109 call this%ksp_init(max_iter)
116 class(
cacg_t),
intent(inout) :: this
120 if (
allocated(this%PR))
then
124 if (
allocated(this%r))
then
128 if (
allocated(this%p))
then
138 function cacg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
139 class(
cacg_t),
intent(inout) :: this
140 class(
ax_t),
intent(inout) :: ax
141 type(
field_t),
intent(inout) :: x
142 integer,
intent(in) :: n
143 real(kind=
rp),
dimension(n),
intent(inout) :: f
144 type(
coef_t),
intent(inout) :: coef
146 type(
gs_t),
intent(inout) :: gs_h
148 integer,
optional,
intent(in) :: niter
149 integer :: i, j, k, l, iter, max_iter, s, ierr, it
150 real(kind=
rp) :: rnorm, rtr, rtz1, tmp
151 real(kind=
rp) :: beta(this%s+1), alpha(this%s+1), alpha1, alpha2, norm_fac
152 real(kind=
rp),
dimension(4*this%s+1,4*this%s+1) :: tt, g, gtt, temp, temp2
153 real(kind=
rp) :: p_c(4*this%s+1,this%s+1)
154 real(kind=
rp) :: r_c(4*this%s+1,this%s+1)
155 real(kind=
rp) :: z_c(4*this%s+1,this%s+1)
156 real(kind=
rp) :: x_c(4*this%s+1,this%s+1)
158 associate(pr => this%PR, r => this%r, p => this%p)
160 if (
present(niter))
then
163 max_iter = this%max_iter
165 norm_fac = 1.0_rp / sqrt(coef%volume)
170 call this%M%solve(p, r, n)
172 rtr =
glsc3(r, coef%mult, r, n)
173 rnorm = sqrt(rtr)*norm_fac
174 ksp_results%res_start = rnorm
175 ksp_results%res_final = rnorm
178 if(
abscmp(rnorm, 0.0_rp))
return
179 call this%monitor_start(
'CACG')
180 do while (iter < max_iter)
183 call copy(pr(1,2*s+2), r, n)
187 if (mod(i,2) .eq. 0)
then
188 call ax%compute(pr(1,i), pr(1,i-1), coef, x%msh, x%Xh)
189 call gs_h%gs_op_vector(pr(1,i), n, gs_op_add)
192 call this%M%solve(pr(1,i), pr(1,i-1), n)
197 if (mod(i,2) == 0)
then
198 call this%M%solve(pr(1,i+1), pr(1,i), n)
200 call ax%compute(pr(1,i+1), pr(1,i), coef, x%msh, x%Xh)
201 call gs_h%gs_op_vector(pr(1,i+1), n, gs_op_add)
207 call rzero(p_c, (4*s+1) * (s+1))
209 call rzero(r_c, (4*s+1) * (s+1))
210 r_c(2*s+2,1) = 1.0_rp
211 call mxm(tt, 4*s+1, r_c, 4*s+1, z_c,s+1)
212 call rzero(x_c, (4*s+1) * (s+1))
213 call rzero(temp, (4*s+1)**2)
215 do i = 0, n, neko_blk_size
217 if (i + neko_blk_size .le. n)
then
221 do k = 1, neko_blk_size
222 temp(it,1) = temp(it,1) &
223 + pr(i+k,j) * pr(i+k,l) * coef%mult(i+k,1,1,1)
232 temp(it,1) = temp(it,1) &
233 + pr(i+k,j) * pr(i+k,l) * coef%mult(i+k,1,1,1)
240 call mpi_allreduce(temp, temp2, it, &
251 call mxm(g,4*s+1, tt, 4*s+1,gtt,4*s+1)
256 call mxm(g, 4*s+1, r_c(1,j), 4*s+1,temp, 1)
257 call mxm(gtt, 4*s+1, p_c(1,j), 4*s+1,temp2, 1)
261 alpha1 = alpha1 + temp(i,1) * z_c(i,j)
262 alpha2 = alpha2 + temp2(i,1) * p_c(i,j)
264 alpha(j) = alpha1/alpha2
267 x_c(i,j+1) = x_c(i,j) + alpha(j) * p_c(i,j)
270 tmp = tmp + tt(i,k) * p_c(k,j)
272 r_c(i,j+1) = r_c(i,j) - alpha(j)*tmp
275 tmp = tmp + tt(i,k)*r_c(k,j+1)
280 call mxm(g,4*s+1,r_c(1,j+1),4*s+1,temp2,1)
283 alpha2 = alpha2 + temp2(i,1)*z_c(i,j+1)
285 beta(j) = alpha2 / alpha1
287 p_c(i,j+1) = z_c(i,j+1) + beta(j)*p_c(i,j)
294 do i = 0, n, neko_blk_size
295 if (i + neko_blk_size .le. n)
then
297 do k = 1, neko_blk_size
298 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + pr(i+k,j) * x_c(j,s+1)
299 p(i+k) = p(i+k) + pr(i+k,j) * p_c(j,s+1)
300 tmp = pr(i+k,j) * r_c(j,s+1)
301 r(i+k) = r(i+k) + tmp
304 do k = 1, neko_blk_size
305 rtr = rtr + r(i+k)**2 * coef%mult(i+k,1,1,1)
310 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + pr(i+k,j) * x_c(j,s+1)
311 p(i+k) = p(i+k) + pr(i+k,j) * p_c(j,s+1)
312 tmp = pr(i+k,j) * r_c(j,s+1)
313 r(i+k) = r(i+k) + tmp
317 rtr = rtr + r(i+k)**2 * coef%mult(i+k,1,1,1)
322 call mpi_allreduce(rtr, tmp, 1, &
324 rnorm = norm_fac*sqrt(tmp)
325 call this%monitor_iter(iter, rnorm)
326 if( rnorm <= this%abs_tol)
exit
328 call this%monitor_stop()
329 ksp_results%res_final = rnorm
330 ksp_results%iter = iter
338 integer,
intent(in) :: s
339 real(kind=rp),
intent(inout) :: tt(4*s+1,4*s+1)
341 mlen = (4*s+1)*(4*s+1)
347 tt(2*s+2+i,2*s+1+i) = 1.0_rp
353 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
354 class(
cacg_t),
intent(inout) :: this
355 class(ax_t),
intent(inout) :: ax
356 type(field_t),
intent(inout) :: x
357 type(field_t),
intent(inout) :: y
358 type(field_t),
intent(inout) :: z
359 integer,
intent(in) :: n
360 real(kind=rp),
dimension(n),
intent(inout) :: fx
361 real(kind=rp),
dimension(n),
intent(inout) :: fy
362 real(kind=rp),
dimension(n),
intent(inout) :: fz
363 type(coef_t),
intent(inout) :: coef
364 type(bc_list_t),
intent(inout) :: blstx
365 type(bc_list_t),
intent(inout) :: blsty
366 type(bc_list_t),
intent(inout) :: blstz
367 type(gs_t),
intent(inout) :: gs_h
368 type(ksp_monitor_t),
dimension(3) :: ksp_results
369 integer,
optional,
intent(in) :: niter
371 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
372 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
373 ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
Defines a Matrix-vector product.
Defines a boundary condition.
subroutine, public bc_list_apply_scalar(bclst, x, n, t, tstep)
Apply a list of boundary conditions to a scalar field.
Defines a communication avoiding Conjugate Gradient method.
type(ksp_monitor_t) function, dimension(3) cacg_solve_coupled(this, Ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
S-step CA PCG coupled solve.
subroutine construct_basis_matrix(Tt, s)
Monomial matrix constuction, not sparse.
subroutine cacg_init(this, n, max_iter, M, s, rel_tol, abs_tol, monitor)
Initialise a s-step CA PCG solver.
type(ksp_monitor_t) function cacg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
S-step CA PCG solve.
subroutine cacg_free(this)
Deallocate a s-step CA PCG solver.
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 x_update(a, b, c, c1, c2, n)
Returns .
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
Wrapper for all matrix-matrix product implementations.
subroutine, public mxm(a, n1, b, n2, c, n3)
Compute matrix-matrix product for contiguously packed matrices A,B, and C.
integer, parameter, public rp
Global precision used in computations.
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Base type for a matrix-vector product providing .
A list of boundary conditions.
S-step communication avoiding 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.