Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
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!
34module 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_list, only : bc_list_t
44 use identity, only : ident_t
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 character(len=18) :: name = ""
60 integer :: iter
62 real(kind=rp) :: res_start
64 real(kind=rp) :: res_final
66 logical :: converged = .false.
67 contains
68 procedure, pass(this) :: print_header => krylov_monitor_print_header
69 procedure, pass(this) :: print_result => krylov_monitor_print_result
70 end type ksp_monitor_t
71
73 type, public, abstract :: ksp_t
74 class(pc_t), pointer :: m => null()
75 real(kind=rp) :: rel_tol
76 real(kind=rp) :: abs_tol
77 integer :: max_iter
78 class(pc_t), allocatable :: m_ident
79 logical :: monitor
80 contains
82 procedure(ksp_init_intrf), deferred, pass(this) :: init
84 procedure, pass(this) :: ksp_init => krylov_init
86 procedure, pass(this) :: ksp_free => krylov_free
88 procedure, pass(this) :: set_pc => krylov_set_pc
90 procedure(ksp_method), pass(this), deferred :: solve
92 procedure(ksp_method_coupled), pass(this), deferred :: solve_coupled
94 procedure, pass(this) :: monitor_start => krylov_monitor_start
96 procedure, pass(this) :: monitor_stop => krylov_monitor_stop
98 procedure, pass(this) :: monitor_iter => krylov_monitor_iter
100 procedure, pass(this) :: is_converged => krylov_is_converged
102 procedure(ksp_t_free), pass(this), deferred :: free
103 end type ksp_t
104
112 abstract interface
113 subroutine ksp_init_intrf(this, n, max_iter, M, rel_tol, abs_tol, monitor)
114 import :: pc_t, ksp_t, rp
115 implicit none
116 class(ksp_t), target, intent(inout) :: this
117 integer, intent(in) :: max_iter
118 class(pc_t), optional, intent(in), target :: M
119 integer, intent(in) :: n
120 real(kind=rp), optional, intent(in) :: rel_tol
121 real(kind=rp), optional, intent(in) :: abs_tol
122 logical, optional, intent(in) :: monitor
123 end subroutine ksp_init_intrf
124 end interface
125
135 abstract interface
136 function ksp_method(this, Ax, x, f, n, coef, blst, gs_h, niter) &
137 result(ksp_results)
138 import :: bc_list_t
139 import :: field_t
140 import :: ksp_t
141 import :: coef_t
142 import :: gs_t
143 import :: ax_t
144 import :: ksp_monitor_t
145 import rp
146 implicit none
147 class(ksp_t), intent(inout) :: this
148 class(ax_t), intent(in) :: ax
149 type(field_t), intent(inout) :: x
150 integer, intent(in) :: n
151 real(kind=rp), dimension(n), intent(in) :: f
152 type(coef_t), intent(inout) :: coef
153 type(bc_list_t), intent(inout) :: blst
154 type(gs_t), intent(inout) :: gs_h
155 integer, optional, intent(in) :: niter
156 type(ksp_monitor_t) :: ksp_results
157 end function ksp_method
158 end interface
159
173 abstract interface
174 function ksp_method_coupled(this, Ax, x, y, z, fx, fy, fz, &
175 n, coef, blstx, blsty, blstz, gs_h, niter) result(ksp_results)
176 import :: bc_list_t
177 import :: field_t
178 import :: ksp_t
179 import :: coef_t
180 import :: gs_t
181 import :: ax_t
182 import :: ksp_monitor_t
183 import rp
184 implicit none
185 class(ksp_t), intent(inout) :: this
186 class(ax_t), intent(in) :: ax
187 type(field_t), intent(inout) :: x
188 type(field_t), intent(inout) :: y
189 type(field_t), intent(inout) :: z
190 integer, intent(in) :: n
191 real(kind=rp), dimension(n), intent(in) :: fx
192 real(kind=rp), dimension(n), intent(in) :: fy
193 real(kind=rp), dimension(n), intent(in) :: fz
194 type(coef_t), intent(inout) :: coef
195 type(bc_list_t), intent(inout) :: blstx
196 type(bc_list_t), intent(inout) :: blsty
197 type(bc_list_t), intent(inout) :: blstz
198 type(gs_t), intent(inout) :: gs_h
199 integer, optional, intent(in) :: niter
200 type(ksp_monitor_t), dimension(3) :: ksp_results
201 end function ksp_method_coupled
202 end interface
203
205 abstract interface
206 subroutine ksp_t_free(this)
207 import :: ksp_t
208 class(ksp_t), intent(inout) :: this
209 end subroutine ksp_t_free
210 end interface
211
212 interface
213
221 module subroutine krylov_solver_factory(object, n, type_name, &
222 max_iter, abstol, m, monitor)
223 class(ksp_t), allocatable, intent(inout) :: object
224 integer, intent(in), value :: n
225 character(len=*), intent(in) :: type_name
226 integer, intent(in) :: max_iter
227 real(kind=rp), optional :: abstol
228 class(pc_t), optional, intent(in), target :: m
229 logical, optional, intent(in) :: monitor
230 end subroutine krylov_solver_factory
231
232 end interface
233
234 public :: krylov_solver_factory
235contains
236
242 subroutine krylov_init(this, max_iter, rel_tol, abs_tol, M, monitor)
243 class(ksp_t), target, intent(inout) :: this
244 integer, intent(in) :: max_iter
245 real(kind=rp), optional, intent(in) :: rel_tol
246 real(kind=rp), optional, intent(in) :: abs_tol
247 class(pc_t), optional, target, intent(in) :: m
248 logical, optional, intent(in) :: monitor
249
250 call krylov_free(this)
251
252 if (present(rel_tol)) then
253 this%rel_tol = rel_tol
254 else
255 this%rel_tol = ksp_rel_tol
256 end if
257
258 if (present(abs_tol)) then
259 this%abs_tol = abs_tol
260 else
261 this%abs_tol = ksp_abs_tol
262 end if
263
264 this%max_iter = max_iter
265
266 if (present(m)) then
267 this%M => m
268 else
269 if (.not. associated(this%M)) then
270 if (neko_bcknd_device .eq. 1) then
271 allocate(device_ident_t::this%M_ident)
272 else
273 allocate(ident_t::this%M_ident)
274 end if
275 this%M => this%M_ident
276 end if
277 end if
278
279 if (present(monitor)) then
280 this%monitor = monitor
281 else
282 this%monitor = .false.
283 end if
284
285 end subroutine krylov_init
286
288 subroutine krylov_free(this)
289 class(ksp_t), intent(inout) :: this
290
292
293 end subroutine krylov_free
294
297 subroutine krylov_set_pc(this, M)
298 class(ksp_t), intent(inout) :: this
299 class(pc_t), target, intent(in) :: M
300
301 if (associated(this%M)) then
302 select type (pc => this%M)
303 type is (ident_t)
304 type is (device_ident_t)
305 class default
306 call neko_error('Preconditioner already defined')
307 end select
308 end if
309
310 this%M => m
311
312 end subroutine krylov_set_pc
313
315 subroutine krylov_monitor_start(this, name)
316 class(ksp_t), intent(in) :: this
317 character(len=*) :: name
318 character(len=LOG_SIZE) :: log_buf
319
320 if (this%monitor) then
321 write(log_buf, '(A)') 'Krylov monitor (' // trim(name) // ')'
322 call neko_log%section(trim(log_buf))
323 call neko_log%newline()
324 call neko_log%begin()
325 write(log_buf, '(A)') ' Iter. Residual'
326 call neko_log%message(log_buf)
327 write(log_buf, '(A)') '-------------------------'
328 call neko_log%message(log_buf)
329 end if
330 end subroutine krylov_monitor_start
331
333 subroutine krylov_monitor_stop(this)
334 class(ksp_t), intent(in) :: this
335
336 if (this%monitor) then
337 call neko_log%end()
338 call neko_log%end_section()
339 call neko_log%newline()
340 end if
341 end subroutine krylov_monitor_stop
342
343
345 subroutine krylov_monitor_iter(this, iter, rnorm)
346 class(ksp_t), intent(in) :: this
347 integer, intent(in) :: iter
348 real(kind=rp), intent(in) :: rnorm
349 character(len=LOG_SIZE) :: log_buf
350
351 if (this%monitor) then
352 write(log_buf, '(I6,E18.9)') iter, rnorm
353 call neko_log%message(log_buf)
354 end if
355
356 end subroutine krylov_monitor_iter
357
366 pure function krylov_is_converged(this, iter, residual) result(converged)
367 class(ksp_t), intent(in) :: this
368 integer, intent(in) :: iter
369 real(kind=rp), intent(in) :: residual
370 logical :: converged
371
372 converged = .true.
373 if (iter .ge. this%max_iter) converged = .false.
374 if (residual .gt. this%abs_tol) converged = .false.
375
376 end function krylov_is_converged
377
379 subroutine krylov_monitor_print_header(this)
380 class(ksp_monitor_t), intent(in) :: this
381 character(len=LOG_SIZE) :: log_buf
382
383 write(log_buf, '((A5,7x),A3,(A5,13x),1x,A6,3x,A15,3x,A15)') &
384 'Step:', ' | ', 'Field:', 'Iters:', &
385 'Start residual:', 'Final residual:'
386 call neko_log%message(log_buf)
387
388 end subroutine krylov_monitor_print_header
389
391 subroutine krylov_monitor_print_result(this, step)
392 class(ksp_monitor_t), intent(in) :: this
393 integer, intent(in) :: step
394 character(len=LOG_SIZE) :: log_buf
395 character(len=12) :: step_str
396 character(len=:), allocatable :: output_format
397
398 if (this%name .eq. "") call neko_error('Krylov solver name is not set')
399
400 ! Define the output format
401 output_format = '(A12,A3,A18,1x,I6,3x,E15.9,3x,E15.9)'
402 write(step_str, '(I12)') step
403 step_str = adjustl(step_str)
404
405 write(log_buf, output_format) &
406 step_str, ' | ' , adjustl(this%name), this%iter, &
407 this%res_start, this%res_final
408 call neko_log%message(log_buf)
409
410 end subroutine krylov_monitor_print_result
411
412end module krylov
__device__ T solve(const T u, const T y, const T guess, const T nu, const T kappa, const T B)
Abstract interface for a Krylov method's constructor.
Definition krylov.f90:113
Abstract interface for a Krylov method's coupled solve routine.
Definition krylov.f90:174
Abstract interface for a Krylov method's solve routine.
Definition krylov.f90:136
Abstract interface for deallocating a Krylov method.
Definition krylov.f90:206
Defines a Matrix-vector product.
Definition ax.f90:34
Defines a list of bc_t.
Definition bc_list.f90:34
Coefficients.
Definition coef.f90:34
Identity Krylov preconditioner for accelerators.
Defines a field.
Definition field.f90:34
Gather-scatter.
Krylov preconditioner (identity)
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:289
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:346
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:243
subroutine krylov_monitor_print_header(this)
Print the Krylov solver's result header.
Definition krylov.f90:380
subroutine krylov_monitor_start(this, name)
Monitor start.
Definition krylov.f90:316
pure logical function krylov_is_converged(this, iter, residual)
Check for convergence.
Definition krylov.f90:367
subroutine krylov_set_pc(this, m)
Setup a Krylov solver's preconditioner.
Definition krylov.f90:298
subroutine krylov_monitor_stop(this)
Monitor stop.
Definition krylov.f90:334
subroutine krylov_monitor_print_result(this, step)
Print the Krylov solver's result.
Definition krylov.f90:392
Logging routines.
Definition log.f90:34
type(log_t), public neko_log
Global log stream.
Definition log.f90:70
integer, parameter, public log_size
Definition log.f90:46
Defines a mesh.
Definition mesh.f90:34
Build configurations.
integer, parameter neko_bcknd_device
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:284
Base type for a matrix-vector product providing .
Definition ax.f90:43
A list of allocatable `bc_t`. Follows the standard interface of lists.
Definition bc_list.f90:48
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.
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:73
Defines a canonical Krylov preconditioner.
Definition precon.f90:40