Neko  0.9.99
A portable framework for high-order spectral element flow simulations
bicgstab.f90
Go to the documentation of this file.
1 ! Copyright (c) 2021-2024, The Neko Authors
2 ! All rights reserved.
3 !
4 ! Redistribution and use in source and binary forms, with or without
5 ! modification, are permitted provided that the following conditions
6 ! are met:
7 !
8 ! * Redistributions of source code must retain the above copyright
9 ! notice, this list of conditions and the following disclaimer.
10 !
11 ! * Redistributions in binary form must reproduce the above
12 ! copyright notice, this list of conditions and the following
13 ! disclaimer in the documentation and/or other materials provided
14 ! with the distribution.
15 !
16 ! * Neither the name of the authors nor the names of its
17 ! contributors may be used to endorse or promote products derived
18 ! from this software without specific prior written permission.
19 !
20 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23 ! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24 ! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25 ! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26 ! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27 ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30 ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 ! POSSIBILITY OF SUCH DAMAGE.
32 !
34 module bicgstab
35  use num_types, only: rp
37  use precon, only : pc_t
38  use ax_product, only : ax_t
39  use field, only : field_t
40  use coefs, only : coef_t
41  use gather_scatter, only : gs_t, gs_op_add
42  use bc, only : bc_list_t, bc_list_apply
43  use math, only : glsc3, rzero, copy, neko_eps, add2s2, x_update, &
45  use utils, only : neko_error
46  implicit none
47  private
48 
50  type, public, extends(ksp_t) :: bicgstab_t
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(:)
58  contains
60  procedure, pass(this) :: init => bicgstab_init
62  procedure, pass(this) :: free => bicgstab_free
64  procedure, pass(this) :: solve => bicgstab_solve
66  procedure, pass(this) :: solve_coupled => bicgstab_solve_coupled
67  end type bicgstab_t
68 
69 contains
70 
72  subroutine bicgstab_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
73  class(bicgstab_t), intent(inout) :: this
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
80 
81 
82  call this%free()
83 
84  allocate(this%p(n))
85  allocate(this%p_hat(n))
86  allocate(this%r(n))
87  allocate(this%s(n))
88  allocate(this%s_hat(n))
89  allocate(this%t(n))
90  allocate(this%v(n))
91  if (present(m)) then
92  this%M => m
93  end if
94 
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)
109  else
110  call this%ksp_init(max_iter)
111  end if
112 
113  end subroutine bicgstab_init
114 
116  subroutine bicgstab_free(this)
117  class(bicgstab_t), intent(inout) :: this
118 
119  call this%ksp_free()
120 
121  if (allocated(this%v)) then
122  deallocate(this%v)
123  end if
124 
125  if (allocated(this%r)) then
126  deallocate(this%r)
127  end if
128 
129  if (allocated(this%t)) then
130  deallocate(this%t)
131  end if
132 
133  if (allocated(this%p)) then
134  deallocate(this%p)
135  end if
136 
137  if (allocated(this%p_hat)) then
138  deallocate(this%p_hat)
139  end if
140 
141  if (allocated(this%s)) then
142  deallocate(this%s)
143  end if
144 
145  if (allocated(this%s_hat)) then
146  deallocate(this%s_hat)
147  end if
148 
149  nullify(this%M)
150 
151 
152  end subroutine bicgstab_free
153 
155  function bicgstab_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_results)
156  class(bicgstab_t), intent(inout) :: this
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
162  type(bc_list_t), intent(in) :: blst
163  type(gs_t), intent(inout) :: gs_h
164  type(ksp_monitor_t) :: ksp_results
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
169 
170  if (present(niter)) then
171  max_iter = niter
172  else
173  max_iter = this%max_iter
174  end if
175  norm_fac = 1.0_rp / sqrt(coef%volume)
176 
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)
179 
180  call rzero(x%x, n)
181  call copy(r, f, n)
182 
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
188  ksp_results%iter = 0
189  if(abscmp(rnorm, 0.0_rp)) return
190  call this%monitor_start('BiCGStab')
191  do iter = 1, max_iter
192 
193  rho_1 = glsc3(r, coef%mult, f ,n)
194 
195  if (abs(rho_1) .lt. neko_eps) then
196  call neko_error('Bi-CGStab rho failure')
197  end if
198 
199  if (iter .eq. 1) then
200  call copy(p, r, n)
201  else
202  beta = (rho_1 / rho_2) * (alpha / omega)
203  call p_update(p, r, v, beta, omega, n)
204  end if
205 
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)
209  call bc_list_apply(blst, v, n)
210  alpha = rho_1 / glsc3(f, coef%mult, v, n)
211  call copy(s, r, 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)
217  exit
218  end if
219 
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)
223  call bc_list_apply(blst, t, n)
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)
227  call copy(r, s, n)
228  call add2s2(r, t, -omega, n)
229 
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
234  exit
235  end if
236 
237  if (omega .lt. neko_eps) then
238  call neko_error('Bi-CGstab omega failure')
239  end if
240  rho_2 = rho_1
241 
242  end do
243  end associate
244  call this%monitor_stop()
245  ksp_results%res_final = rnorm
246  ksp_results%iter = iter
247  end function bicgstab_solve
248 
250  function bicgstab_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
251  n, coef, blstx, blsty, blstz, gs_h, niter) result(ksp_results)
252  class(bicgstab_t), intent(inout) :: this
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
268 
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)
272 
273  end function bicgstab_solve_coupled
274 
275 end module bicgstab
276 
277 
Defines a Matrix-vector product.
Definition: ax.f90:34
Defines a boundary condition.
Definition: bc.f90:34
Defines various Bi-Conjugate Gradient Stabilized methods.
Definition: bicgstab.f90:34
subroutine bicgstab_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
Constructor.
Definition: bicgstab.f90:73
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.
Definition: bicgstab.f90:252
subroutine bicgstab_free(this)
Deallocate a standard BiCGSTAB solver.
Definition: bicgstab.f90:117
type(ksp_monitor_t) function bicgstab_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Bi-Conjugate Gradient Stabilized method solve.
Definition: bicgstab.f90:156
Coefficients.
Definition: coef.f90:34
Defines a field.
Definition: field.f90:34
Gather-scatter.
Implements the base abstract type for Krylov solvers plus helper types.
Definition: krylov.f90:34
integer, parameter, public ksp_max_iter
Maximum number of iters.
Definition: krylov.f90:51
Definition: math.f90:60
real(kind=rp) function, public glsc3(a, b, c, n)
Weighted inner product .
Definition: math.f90:895
subroutine, public x_update(a, b, c, c1, c2, n)
Returns .
Definition: math.f90:861
subroutine, public copy(a, b, n)
Copy a vector .
Definition: math.f90:239
real(kind=rp), parameter, public neko_eps
Machine epsilon .
Definition: math.f90:70
subroutine, public rzero(a, n)
Zero a real vector.
Definition: math.f90:195
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
Definition: math.f90:673
subroutine, public p_update(a, b, c, c1, c2, n)
Returns .
Definition: math.f90:846
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Krylov preconditioner.
Definition: precon.f90:34
Utilities.
Definition: utils.f90:35
Base type for a matrix-vector product providing .
Definition: ax.f90:43
A list of boundary conditions.
Definition: bc.f90:104
Standard preconditioned Bi-Conjugate Gradient Stabilized method.
Definition: bicgstab.f90:50
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
Type for storing initial and final residuals in a Krylov solver.
Definition: krylov.f90:56
Base abstract type for a canonical Krylov method, solving .
Definition: krylov.f90:66
Defines a canonical Krylov preconditioner.
Definition: precon.f90:40