Neko 0.9.99
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
gmres.f90
Go to the documentation of this file.
1! Copyright (c) 2020-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!
34module gmres
35 use krylov, only : ksp_t, ksp_monitor_t
36 use precon, only : pc_t
37 use ax_product, only : ax_t
38 use num_types, only: rp, xp
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_list, only : bc_list_t
43 use math, only : glsc3, rzero, rone, copy, sub2, cmult2, abscmp
44 use neko_config, only : neko_blk_size
45 use comm
46 implicit none
47 private
48
50 type, public, extends(ksp_t) :: gmres_t
51 integer :: lgmres
52 real(kind=rp), allocatable :: w(:)
53 real(kind=rp), allocatable :: r(:)
54 real(kind=rp), allocatable :: z(:,:)
55 real(kind=rp), allocatable :: v(:,:)
57 real(kind=xp), allocatable :: h(:,:)
58 real(kind=xp), allocatable :: s(:)
59 real(kind=xp), allocatable :: gam(:)
60 real(kind=xp), allocatable :: c(:)
61 contains
62 procedure, pass(this) :: init => gmres_init
63 procedure, pass(this) :: free => gmres_free
64 procedure, pass(this) :: solve => gmres_solve
65 procedure, pass(this) :: solve_coupled => gmres_solve_coupled
66 end type gmres_t
67
68contains
69
71 subroutine gmres_init(this, n, max_iter, M, lgmres, &
72 rel_tol, abs_tol, monitor)
73 class(gmres_t), intent(inout) :: this
74 integer, intent(in) :: n
75 integer, intent(in) :: max_iter
76 class(pc_t), optional, intent(in), target :: M
77 integer, optional, intent(in) :: lgmres
78 real(kind=rp), optional, intent(in) :: rel_tol
79 real(kind=rp), optional, intent(in) :: abs_tol
80 logical, optional, intent(in) :: monitor
81
82 if (present(lgmres)) then
83 this%lgmres = lgmres
84 else
85 this%lgmres = 30
86 end if
87
88
89 call this%free()
90
91 if (present(m)) then
92 this%M => m
93 end if
94
95 allocate(this%w(n))
96 allocate(this%r(n))
97
98 allocate(this%c(this%lgmres))
99 allocate(this%s(this%lgmres))
100 allocate(this%gam(this%lgmres + 1))
101
102 allocate(this%z(n, this%lgmres))
103 allocate(this%v(n, this%lgmres))
104
105 allocate(this%h(this%lgmres, this%lgmres))
106
107 if (present(rel_tol) .and. present(abs_tol) .and. present(monitor)) then
108 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
109 else if (present(rel_tol) .and. present(abs_tol)) then
110 call this%ksp_init(max_iter, rel_tol, abs_tol)
111 else if (present(monitor) .and. present(abs_tol)) then
112 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
113 else if (present(rel_tol) .and. present(monitor)) then
114 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
115 else if (present(rel_tol)) then
116 call this%ksp_init(max_iter, rel_tol = rel_tol)
117 else if (present(abs_tol)) then
118 call this%ksp_init(max_iter, abs_tol = abs_tol)
119 else if (present(monitor)) then
120 call this%ksp_init(max_iter, monitor = monitor)
121 else
122 call this%ksp_init(max_iter)
123 end if
124
125 end subroutine gmres_init
126
128 subroutine gmres_free(this)
129 class(gmres_t), intent(inout) :: this
130
131 call this%ksp_free()
132
133 if (allocated(this%w)) then
134 deallocate(this%w)
135 end if
136
137 if (allocated(this%c)) then
138 deallocate(this%c)
139 end if
140
141 if (allocated(this%r)) then
142 deallocate(this%r)
143 end if
144
145 if (allocated(this%z)) then
146 deallocate(this%z)
147 end if
148
149 if (allocated(this%h)) then
150 deallocate(this%h)
151 end if
152
153 if (allocated(this%v)) then
154 deallocate(this%v)
155 end if
156
157 if (allocated(this%s)) then
158 deallocate(this%s)
159 end if
160
161
162 if (allocated(this%gam)) then
163 deallocate(this%gam)
164 end if
165
166 nullify(this%M)
167
168 end subroutine gmres_free
169
171 function gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
172 result(ksp_results)
173 class(gmres_t), intent(inout) :: this
174 class(ax_t), intent(in) :: ax
175 type(field_t), intent(inout) :: x
176 integer, intent(in) :: n
177 real(kind=rp), dimension(n), intent(in) :: f
178 type(coef_t), intent(inout) :: coef
179 type(bc_list_t), intent(inout) :: blst
180 type(gs_t), intent(inout) :: gs_h
181 type(ksp_monitor_t) :: ksp_results
182 integer, optional, intent(in) :: niter
183 integer :: iter, max_iter
184 integer :: i, j, k, l, ierr
185 real(kind=xp) :: w_plus(neko_blk_size), x_plus(neko_blk_size)
186 real(kind=xp) :: alpha, lr, alpha2, norm_fac
187 real(kind=rp) :: temp, rnorm
188 logical :: conv
189
190 conv = .false.
191 iter = 0
192
193 if (present(niter)) then
194 max_iter = niter
195 else
196 max_iter = this%max_iter
197 end if
198
199 associate(w => this%w, c => this%c, r => this%r, z => this%z, h => this%h, &
200 v => this%v, s => this%s, gam => this%gam)
201
202 norm_fac = 1.0_rp / sqrt(coef%volume)
203 call rzero(x%x, n)
204 gam = 0.0_xp
205 s = 1.0_xp
206 c = 1.0_xp
207 h = 0.0_xp
208 call this%monitor_start('GMRES')
209 do while (.not. conv .and. iter .lt. max_iter)
210
211 if (iter .eq. 0) then
212 call copy(r, f, n)
213 else
214 call copy(r, f, n)
215 call ax%compute(w, x%x, coef, x%msh, x%Xh)
216 call gs_h%op(w, n, gs_op_add)
217 call blst%apply(w, n)
218 call sub2(r, w, n)
219 end if
220
221 gam(1) = sqrt(glsc3(r, r, coef%mult, n))
222 if (iter .eq. 0) then
223 ksp_results%res_start = gam(1) * norm_fac
224 end if
225
226 if (abscmp(gam(1), 0.0_xp)) return
227
228 rnorm = 0.0_rp
229 temp = 1.0_rp / gam(1)
230 call cmult2(v(1,1), r, temp, n)
231 do j = 1, this%lgmres
232 iter = iter+1
233
234 call this%M%solve(z(1,j), v(1,j), n)
235
236 call ax%compute(w, z(1,j), coef, x%msh, x%Xh)
237 call gs_h%op(w, n, gs_op_add)
238 call blst%apply(w, n)
239
240 do l = 1, j
241 h(l,j) = 0.0_rp
242 end do
243
244 do i = 0, n, neko_blk_size
245 if (i + neko_blk_size .le. n) then
246 do l = 1, j
247 do k = 1, neko_blk_size
248 h(l,j) = h(l,j) + &
249 w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
250 end do
251 end do
252 else
253 do k = 1, n-i
254 do l = 1, j
255 h(l,j) = h(l,j) + &
256 w(i+k) * v(i+k,l) * coef%mult(i+k,1,1,1)
257 end do
258 end do
259 end if
260 end do
261
262 call mpi_allreduce(mpi_in_place, h(1,j), j, &
263 mpi_extra_precision, mpi_sum, neko_comm, ierr)
264
265 alpha2 = 0.0_rp
266 do i = 0, n, neko_blk_size
267 if (i + neko_blk_size .le. n) then
268 do k = 1, neko_blk_size
269 w_plus(k) = 0.0_rp
270 end do
271 do l = 1,j
272 do k = 1, neko_blk_size
273 w_plus(k) = w_plus(k) - h(l,j) * v(i+k,l)
274 end do
275 end do
276 do k = 1, neko_blk_size
277 w(i+k) = w(i+k) + w_plus(k)
278 alpha2 = alpha2 + w(i+k)**2 * coef%mult(i+k,1,1,1)
279 end do
280 else
281 do k = 1, n-i
282 w_plus(1) = 0.0_rp
283 do l = 1, j
284 w_plus(1) = w_plus(1) - h(l,j) * v(i+k,l)
285 end do
286 w(i+k) = w(i+k) + w_plus(1)
287 alpha2 = alpha2 + (w(i+k)**2) * coef%mult(i+k,1,1,1)
288 end do
289 end if
290 end do
291
292 call mpi_allreduce(mpi_in_place,alpha2, 1, &
293 mpi_extra_precision, mpi_sum, neko_comm, ierr)
294 alpha = sqrt(alpha2)
295 do i = 1, j-1
296 temp = h(i,j)
297 h(i,j) = c(i)*temp + s(i) * h(i+1,j)
298 h(i+1,j) = -s(i)*temp + c(i) * h(i+1,j)
299 end do
300
301 rnorm = 0.0_rp
302 if (abscmp(alpha, 0.0_xp)) then
303 conv = .true.
304 exit
305 end if
306
307 lr = sqrt(h(j,j) * h(j,j) + alpha2)
308 temp = 1.0_rp / lr
309 c(j) = h(j,j) * temp
310 s(j) = alpha * temp
311 h(j,j) = lr
312 gam(j+1) = -s(j) * gam(j)
313 gam(j) = c(j) * gam(j)
314 rnorm = abs(gam(j+1)) * norm_fac
315 call this%monitor_iter(iter, rnorm)
316 if (rnorm .lt. this%abs_tol) then
317 conv = .true.
318 exit
319 end if
320
321 if (iter + 1 .gt. max_iter) exit
322
323 if (j .lt. this%lgmres) then
324 temp = 1.0_rp / alpha
325 call cmult2(v(1,j+1), w, temp, n)
326 end if
327
328 end do
329
330 j = min(j, this%lgmres)
331 do k = j, 1, -1
332 temp = gam(k)
333 do i = j, k+1, -1
334 temp = temp - h(k,i) * c(i)
335 end do
336 c(k) = temp / h(k,k)
337 end do
338
339 do i = 0, n, neko_blk_size
340 if (i + neko_blk_size .le. n) then
341 do k = 1, neko_blk_size
342 x_plus(k) = 0.0_rp
343 end do
344 do l = 1,j
345 do k = 1, neko_blk_size
346 x_plus(k) = x_plus(k) + c(l) * z(i+k,l)
347 end do
348 end do
349 do k = 1, neko_blk_size
350 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(k)
351 end do
352 else
353 do k = 1, n-i
354 x_plus(1) = 0.0_rp
355 do l = 1, j
356 x_plus(1) = x_plus(1) + c(l) * z(i+k,l)
357 end do
358 x%x(i+k,1,1,1) = x%x(i+k,1,1,1) + x_plus(1)
359 end do
360 end if
361 end do
362 end do
363
364 end associate
365 call this%monitor_stop()
366 ksp_results%res_final = rnorm
367 ksp_results%iter = iter
368 ksp_results%converged = this%is_converged(iter, rnorm)
369
370 end function gmres_solve
371
373 function gmres_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
374 n, coef, blstx, blsty, blstz, gs_h, niter) result(ksp_results)
375 class(gmres_t), intent(inout) :: this
376 class(ax_t), intent(in) :: ax
377 type(field_t), intent(inout) :: x
378 type(field_t), intent(inout) :: y
379 type(field_t), intent(inout) :: z
380 integer, intent(in) :: n
381 real(kind=rp), dimension(n), intent(in) :: fx
382 real(kind=rp), dimension(n), intent(in) :: fy
383 real(kind=rp), dimension(n), intent(in) :: fz
384 type(coef_t), intent(inout) :: coef
385 type(bc_list_t), intent(inout) :: blstx
386 type(bc_list_t), intent(inout) :: blsty
387 type(bc_list_t), intent(inout) :: blstz
388 type(gs_t), intent(inout) :: gs_h
389 type(ksp_monitor_t), dimension(3) :: ksp_results
390 integer, optional, intent(in) :: niter
391
392 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
393 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
394 ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
395
396 end function gmres_solve_coupled
397
398end module gmres
399
400
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
Definition comm.F90:1
type(mpi_comm) neko_comm
MPI communicator.
Definition comm.F90:16
type(mpi_datatype) mpi_extra_precision
Definition comm.F90:24
Defines a field.
Definition field.f90:34
Gather-scatter.
Defines various GMRES methods.
Definition gmres.f90:34
subroutine gmres_init(this, n, max_iter, m, lgmres, rel_tol, abs_tol, monitor)
Initialise a standard GMRES solver.
Definition gmres.f90:73
type(ksp_monitor_t) function gmres_solve(this, ax, x, f, n, coef, blst, gs_h, niter)
Standard GMRES solve.
Definition gmres.f90:173
type(ksp_monitor_t) function, dimension(3) gmres_solve_coupled(this, ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
Standard GMRES coupled solve.
Definition gmres.f90:375
subroutine gmres_free(this)
Deallocate a standard GMRES solver.
Definition gmres.f90:129
Implements the base abstract type for Krylov solvers plus helper types.
Definition krylov.f90:34
Definition math.f90:60
subroutine, public cmult2(a, b, c, n)
Multiplication by constant c .
Definition math.f90:700
real(kind=rp) function, public glsc3(a, b, c, n)
Weighted inner product .
Definition math.f90:894
subroutine, public rone(a, n)
Set all elements to one.
Definition math.f90:227
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:238
subroutine, public rzero(a, n)
Zero a real vector.
Definition math.f90:194
subroutine, public sub2(a, b, n)
Vector substraction .
Definition math.f90:628
Build configurations.
integer, parameter neko_blk_size
integer, parameter, public xp
Definition num_types.f90:14
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Krylov preconditioner.
Definition precon.f90:34
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:47
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
Standard preconditioned generalized minimal residual method.
Definition gmres.f90:50
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:68
Defines a canonical Krylov preconditioner.
Definition precon.f90:40