51 real(kind=
rp),
allocatable :: p(:)
52 real(kind=
rp),
allocatable :: p_hat(:)
53 real(kind=
rp),
allocatable :: r(:)
54 real(kind=
rp),
allocatable :: s(:)
55 real(kind=
rp),
allocatable :: s_hat(:)
56 real(kind=
rp),
allocatable :: t(:)
57 real(kind=
rp),
allocatable :: v(:)
72 subroutine bicgstab_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
74 class(
pc_t),
optional,
intent(in),
target :: M
75 integer,
intent(in) :: n
76 integer,
intent(in) :: max_iter
77 real(kind=
rp),
optional,
intent(in) :: rel_tol
78 real(kind=
rp),
optional,
intent(in) :: abs_tol
79 logical,
optional,
intent(in) :: monitor
85 allocate(this%p_hat(n))
88 allocate(this%s_hat(n))
95 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
96 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
97 else if (
present(rel_tol) .and.
present(abs_tol))
then
98 call this%ksp_init(max_iter, rel_tol, abs_tol)
99 else if (
present(monitor) .and.
present(abs_tol))
then
100 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
101 else if (
present(rel_tol) .and.
present(monitor))
then
102 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
103 else if (
present(rel_tol))
then
104 call this%ksp_init(max_iter, rel_tol = rel_tol)
105 else if (
present(abs_tol))
then
106 call this%ksp_init(max_iter, abs_tol = abs_tol)
107 else if (
present(monitor))
then
108 call this%ksp_init(max_iter, monitor = monitor)
110 call this%ksp_init(max_iter)
121 if (
allocated(this%v))
then
125 if (
allocated(this%r))
then
129 if (
allocated(this%t))
then
133 if (
allocated(this%p))
then
137 if (
allocated(this%p_hat))
then
138 deallocate(this%p_hat)
141 if (
allocated(this%s))
then
145 if (
allocated(this%s_hat))
then
146 deallocate(this%s_hat)
155 function bicgstab_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
result(ksp_results)
157 class(
ax_t),
intent(in) :: ax
158 type(
field_t),
intent(inout) :: x
159 integer,
intent(in) :: n
160 real(kind=
rp),
dimension(n),
intent(in) :: f
161 type(
coef_t),
intent(inout) :: coef
163 type(
gs_t),
intent(inout) :: gs_h
165 integer,
optional,
intent(in) :: niter
166 integer :: iter, max_iter
167 real(kind=
rp) :: rnorm, rtr, norm_fac, gamma
168 real(kind=
rp) :: beta, alpha, omega, rho_1, rho_2
170 if (
present(niter))
then
173 max_iter = this%max_iter
175 norm_fac = 1.0_rp / sqrt(coef%volume)
177 associate(r => this%r, t => this%t, s => this%s, v => this%v, p => this%p, &
178 s_hat => this%s_hat, p_hat => this%p_hat)
183 rtr = sqrt(
glsc3(r, coef%mult, r, n))
184 rnorm = rtr * norm_fac
185 gamma = rnorm * this%rel_tol
186 ksp_results%res_start = rnorm
187 ksp_results%res_final = rnorm
189 if(
abscmp(rnorm, 0.0_rp))
return
190 call this%monitor_start(
'BiCGStab')
191 do iter = 1, max_iter
193 rho_1 =
glsc3(r, coef%mult, f ,n)
199 if (iter .eq. 1)
then
202 beta = (rho_1 / rho_2) * (alpha / omega)
203 call p_update(p, r, v, beta, omega, n)
206 call this%M%solve(p_hat, p, n)
207 call ax%compute(v, p_hat, coef, x%msh, x%Xh)
208 call gs_h%op(v, n, gs_op_add)
210 alpha = rho_1 /
glsc3(f, coef%mult, v, n)
212 call add2s2(s, v, -alpha, n)
213 rtr =
glsc3(s, coef%mult, s, n)
214 rnorm = sqrt(rtr) * norm_fac
215 if (rnorm .lt. this%abs_tol .or. rnorm .lt. gamma)
then
216 call add2s2(x%x, p_hat, alpha,n)
220 call this%M%solve(s_hat, s, n)
221 call ax%compute(t, s_hat, coef, x%msh, x%Xh)
222 call gs_h%op(t, n, gs_op_add)
224 omega =
glsc3(t, coef%mult, s, n) &
225 /
glsc3(t, coef%mult, t, n)
226 call x_update(x%x, p_hat, s_hat, alpha, omega, n)
228 call add2s2(r, t, -omega, n)
230 rtr =
glsc3(r, coef%mult, r, n)
231 rnorm = sqrt(rtr) * norm_fac
232 call this%monitor_iter(iter, rnorm)
233 if (rnorm .lt. this%abs_tol .or. rnorm .lt. gamma)
then
244 call this%monitor_stop()
245 ksp_results%res_final = rnorm
246 ksp_results%iter = iter
251 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
253 class(ax_t),
intent(in) :: ax
254 type(field_t),
intent(inout) :: x
255 type(field_t),
intent(inout) :: y
256 type(field_t),
intent(inout) :: z
257 integer,
intent(in) :: n
258 real(kind=rp),
dimension(n),
intent(in) :: fx
259 real(kind=rp),
dimension(n),
intent(in) :: fy
260 real(kind=rp),
dimension(n),
intent(in) :: fz
261 type(coef_t),
intent(inout) :: coef
262 type(bc_list_t),
intent(in) :: blstx
263 type(bc_list_t),
intent(in) :: blsty
264 type(bc_list_t),
intent(in) :: blstz
265 type(gs_t),
intent(inout) :: gs_h
266 type(ksp_monitor_t),
dimension(3) :: ksp_results
267 integer,
optional,
intent(in) :: niter
269 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
270 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
271 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 Bi-Conjugate Gradient Stabilized methods.
subroutine bicgstab_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
Constructor.
type(ksp_monitor_t) function, dimension(3) bicgstab_solve_coupled(this, Ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
Standard BiCGSTAB coupled solve.
subroutine bicgstab_free(this)
Deallocate a standard BiCGSTAB solver.
type(ksp_monitor_t) function bicgstab_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Bi-Conjugate Gradient Stabilized method solve.
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 .
real(kind=rp), parameter, public neko_eps
Machine epsilon .
subroutine, public rzero(a, n)
Zero a real vector.
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
subroutine, public p_update(a, b, c, c1, c2, n)
Returns .
integer, parameter, public rp
Global precision used in computations.
Base type for a matrix-vector product providing .
A list of boundary conditions.
Standard preconditioned Bi-Conjugate Gradient Stabilized 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.