Neko  0.9.0
A portable framework for high-order spectral element flow simulations
krylov.f90
Go to the documentation of this file.
1 ! Copyright (c) 2020-2023, 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 krylov
35  use gather_scatter, only : gs_t, gs_op_add
36  use ax_product, only : ax_t
37  use num_types, only: rp, c_rp
38  use precon, only : pc_t
39  use coefs, only : coef_t
40  use mesh, only : mesh_t
41  use field, only : field_t
42  use utils, only : neko_error, neko_warning
43  use bc, only : bc_list_t
44  use identity, only : ident_t
46  use neko_config, only : neko_bcknd_device
47  use logger, only : neko_log, log_size
48  implicit none
49  private
50 
51  integer, public, parameter :: ksp_max_iter = 1e3
52  real(kind=rp), public, parameter :: ksp_abs_tol = 1d-9
53  real(kind=rp), public, parameter :: ksp_rel_tol = 1d-9
54 
56  type, public :: ksp_monitor_t
58  integer :: iter
60  real(kind=rp) :: res_start
62  real(kind=rp) :: res_final
63  end type ksp_monitor_t
64 
66  type, public, abstract :: ksp_t
67  class(pc_t), pointer :: m => null()
68  real(kind=rp) :: rel_tol
69  real(kind=rp) :: abs_tol
70  integer :: max_iter
71  class(pc_t), allocatable :: m_ident
72  logical :: monitor
73  contains
75  procedure, pass(this) :: ksp_init => krylov_init
77  procedure, pass(this) :: ksp_free => krylov_free
79  procedure, pass(this) :: set_pc => krylov_set_pc
81  procedure(ksp_method), pass(this), deferred :: solve
83  procedure(ksp_method_coupled), pass(this), deferred :: solve_coupled
85  procedure, pass(this) :: monitor_start => krylov_monitor_start
87  procedure, pass(this) :: monitor_stop => krylov_monitor_stop
89  procedure, pass(this) :: monitor_iter => krylov_monitor_iter
91  procedure(ksp_t_free), pass(this), deferred :: free
92  end type ksp_t
93 
94 
104  abstract interface
105  function ksp_method(this, Ax, x, f, n, coef, blst, gs_h, niter) &
106  result(ksp_results)
107  import :: bc_list_t
108  import :: field_t
109  import :: ksp_t
110  import :: coef_t
111  import :: gs_t
112  import :: ax_t
113  import :: ksp_monitor_t
114  import rp
115  implicit none
116  class(ksp_t), intent(inout) :: this
117  class(ax_t), intent(inout) :: ax
118  type(field_t), intent(inout) :: x
119  integer, intent(in) :: n
120  real(kind=rp), dimension(n), intent(inout) :: f
121  type(coef_t), intent(inout) :: coef
122  type(bc_list_t), intent(inout) :: blst
123  type(gs_t), intent(inout) :: gs_h
124  integer, optional, intent(in) :: niter
125  type(ksp_monitor_t) :: ksp_results
126  end function ksp_method
127  end interface
128 
142  abstract interface
143  function ksp_method_coupled(this, Ax, x, y, z, fx, fy, fz, &
144  n, coef, blstx, blsty, blstz, gs_h, niter) result(ksp_results)
145  import :: bc_list_t
146  import :: field_t
147  import :: ksp_t
148  import :: coef_t
149  import :: gs_t
150  import :: ax_t
151  import :: ksp_monitor_t
152  import rp
153  implicit none
154  class(ksp_t), intent(inout) :: this
155  class(ax_t), intent(inout) :: ax
156  type(field_t), intent(inout) :: x
157  type(field_t), intent(inout) :: y
158  type(field_t), intent(inout) :: z
159  integer, intent(in) :: n
160  real(kind=rp), dimension(n), intent(inout) :: fx
161  real(kind=rp), dimension(n), intent(inout) :: fy
162  real(kind=rp), dimension(n), intent(inout) :: fz
163  type(coef_t), intent(inout) :: coef
164  type(bc_list_t), intent(inout) :: blstx
165  type(bc_list_t), intent(inout) :: blsty
166  type(bc_list_t), intent(inout) :: blstz
167  type(gs_t), intent(inout) :: gs_h
168  integer, optional, intent(in) :: niter
169  type(ksp_monitor_t), dimension(3) :: ksp_results
170  end function ksp_method_coupled
171  end interface
172 
174  abstract interface
175  subroutine ksp_t_free(this)
176  import :: ksp_t
177  class(ksp_t), intent(inout) :: this
178  end subroutine ksp_t_free
179  end interface
180 
181  interface
182 
190  module subroutine krylov_solver_factory(object, n, type_name, &
191  max_iter, abstol, m, monitor)
192  class(ksp_t), allocatable, target, intent(inout) :: object
193  integer, intent(in), value :: n
194  character(len=*), intent(in) :: type_name
195  integer, intent(in) :: max_iter
196  real(kind=rp), optional :: abstol
197  class(pc_t), optional, intent(inout), target :: m
198  logical, optional, intent(in) :: monitor
199  end subroutine krylov_solver_factory
200 
202  module subroutine krylov_solver_destroy(object)
203  class(ksp_t), allocatable, intent(inout) :: object
204  end subroutine krylov_solver_destroy
205  end interface
206 
207  public :: krylov_solver_factory, krylov_solver_destroy
208 contains
209 
215  subroutine krylov_init(this, max_iter, rel_tol, abs_tol, M, monitor)
216  class(ksp_t), target, intent(inout) :: this
217  integer, intent(in) :: max_iter
218  real(kind=rp), optional, intent(in) :: rel_tol
219  real(kind=rp), optional, intent(in) :: abs_tol
220  class(pc_t), optional, target, intent(in) :: m
221  logical, optional, intent(in) :: monitor
222 
223  call krylov_free(this)
224 
225  if (present(rel_tol)) then
226  this%rel_tol = rel_tol
227  else
228  this%rel_tol = ksp_rel_tol
229  end if
230 
231  if (present(abs_tol)) then
232  this%abs_tol = abs_tol
233  else
234  this%abs_tol = ksp_abs_tol
235  end if
236 
237  this%max_iter = max_iter
238 
239  if (present(m)) then
240  this%M => m
241  else
242  if (.not. associated(this%M)) then
243  if (neko_bcknd_device .eq. 1) then
244  allocate(device_ident_t::this%M_ident)
245  else
246  allocate(ident_t::this%M_ident)
247  end if
248  this%M => this%M_ident
249  end if
250  end if
251 
252  if (present(monitor)) then
253  this%monitor = monitor
254  else
255  this%monitor = .false.
256  end if
257 
258  end subroutine krylov_init
259 
261  subroutine krylov_free(this)
262  class(ksp_t), intent(inout) :: this
263 
265 
266  end subroutine krylov_free
267 
270  subroutine krylov_set_pc(this, M)
271  class(ksp_t), intent(inout) :: this
272  class(pc_t), target, intent(in) :: M
273 
274  if (associated(this%M)) then
275  select type (pc => this%M)
276  type is (ident_t)
277  type is (device_ident_t)
278  class default
279  call neko_error('Preconditioner already defined')
280  end select
281  end if
282 
283  this%M => m
284 
285  end subroutine krylov_set_pc
286 
288  subroutine krylov_monitor_start(this, name)
289  class(ksp_t), intent(in) :: this
290  character(len=*) :: name
291  character(len=LOG_SIZE) :: log_buf
292 
293  if (this%monitor) then
294  write(log_buf, '(A)') 'Krylov monitor (' // trim(name) // ')'
295  call neko_log%section(trim(log_buf))
296  call neko_log%newline()
297  call neko_log%begin()
298  write(log_buf, '(A)') ' Iter. Residual'
299  call neko_log%message(log_buf)
300  write(log_buf, '(A)') '---------------------'
301  call neko_log%message(log_buf)
302  end if
303  end subroutine krylov_monitor_start
304 
306  subroutine krylov_monitor_stop(this)
307  class(ksp_t), intent(in) :: this
308 
309  if (this%monitor) then
310  call neko_log%end()
311  call neko_log%end_section()
312  call neko_log%newline()
313  end if
314  end subroutine krylov_monitor_stop
315 
316 
318  subroutine krylov_monitor_iter(this, iter, rnorm)
319  class(ksp_t), intent(in) :: this
320  integer, intent(in) :: iter
321  real(kind=rp), intent(in) :: rnorm
322  character(len=LOG_SIZE) :: log_buf
323 
324  if (this%monitor) then
325  write(log_buf, '(I6,E15.7)') iter, rnorm
326  call neko_log%message(log_buf)
327  end if
328 
329  end subroutine krylov_monitor_iter
330 
331 end module krylov
Abstract interface for a Krylov method's coupled solve routine.
Definition: krylov.f90:143
Abstract interface for a Krylov method's solve routine.
Definition: krylov.f90:105
Abstract interface for deallocating a Krylov method.
Definition: krylov.f90:175
Defines a Matrix-vector product.
Definition: ax.f90:34
Defines a boundary condition.
Definition: bc.f90:34
Coefficients.
Definition: coef.f90:34
Identity Krylov preconditioner for accelerators.
Defines a field.
Definition: field.f90:34
Gather-scatter.
Krylov preconditioner (identity)
Definition: pc_identity.f90:34
Implements the base abstract type for Krylov solvers plus helper types.
Definition: krylov.f90:34
real(kind=rp), parameter, public ksp_rel_tol
Relative tolerance.
Definition: krylov.f90:53
real(kind=rp), parameter, public ksp_abs_tol
Absolut tolerance.
Definition: krylov.f90:52
subroutine krylov_free(this)
Deallocate a Krylov solver.
Definition: krylov.f90:262
integer, parameter, public ksp_max_iter
Maximum number of iters.
Definition: krylov.f90:51
subroutine krylov_monitor_iter(this, iter, rnorm)
Monitor iteration.
Definition: krylov.f90:319
subroutine krylov_monitor_start(this, name)
Monitor start.
Definition: krylov.f90:289
subroutine krylov_init(this, max_iter, rel_tol, abs_tol, M, monitor)
Factory for Krylov solvers. Both creates and initializes the object.
Definition: krylov.f90:216
subroutine krylov_set_pc(this, M)
Setup a Krylov solver's preconditioner.
Definition: krylov.f90:271
subroutine krylov_monitor_stop(this)
Monitor stop.
Definition: krylov.f90:307
Logging routines.
Definition: log.f90:34
type(log_t), public neko_log
Global log stream.
Definition: log.f90:65
integer, parameter, public log_size
Definition: log.f90:42
Defines a mesh.
Definition: mesh.f90:34
Build configurations.
Definition: neko_config.f90:34
integer, parameter neko_bcknd_device
Definition: neko_config.f90:44
integer, parameter, public c_rp
Definition: num_types.f90:13
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
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Definition: utils.f90:245
Base type for a matrix-vector product providing .
Definition: ax.f90:43
A list of boundary conditions.
Definition: bc.f90:104
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
Defines a canonical Krylov preconditioner for accelerators.
Defines a canonical Krylov preconditioner.
Definition: pc_identity.f90:42
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