Neko  0.9.99
A portable framework for high-order spectral element flow simulations
gmres_sx.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 gmres_sx
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
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, rone, copy, cmult2, col2, col3, add2s2, abscmp
44  use comm
45  implicit none
46  private
47 
49  type, public, extends(ksp_t) :: sx_gmres_t
50  integer :: lgmres
51  real(kind=rp), allocatable :: w(:)
52  real(kind=rp), allocatable :: c(:)
53  real(kind=rp), allocatable :: r(:)
54  real(kind=rp), allocatable :: z(:,:)
55  real(kind=rp), allocatable :: h(:,:)
56  real(kind=rp), allocatable :: ml(:)
57  real(kind=rp), allocatable :: v(:,:)
58  real(kind=rp), allocatable :: s(:)
59  real(kind=rp), allocatable :: mu(:)
60  real(kind=rp), allocatable :: gam(:)
61  real(kind=rp), allocatable :: wk1(:)
62  real(kind=rp) :: rnorm
63  contains
64  procedure, pass(this) :: init => sx_gmres_init
65  procedure, pass(this) :: free => sx_gmres_free
66  procedure, pass(this) :: solve => sx_gmres_solve
67  procedure, pass(this) :: solve_coupled => sx_gmres_solve_coupled
68  end type sx_gmres_t
69 
70 contains
71 
73  subroutine sx_gmres_init(this, n, max_iter, M, lgmres, &
74  rel_tol, abs_tol, monitor)
75  class(sx_gmres_t), intent(inout) :: this
76  integer, intent(in) :: n
77  integer, intent(in) :: max_iter
78  class(pc_t), optional, intent(in), target :: M
79  integer, optional, intent(in) :: lgmres
80  real(kind=rp), optional, intent(in) :: rel_tol
81  real(kind=rp), optional, intent(in) :: abs_tol
82  logical, optional, intent(in) :: monitor
83 
84  if (present(lgmres)) then
85  this%lgmres = lgmres
86  else
87  this%lgmres = 30
88  end if
89 
90 
91  call this%free()
92 
93  if (present(m)) then
94  this%M => m
95  end if
96 
97  allocate(this%w(n))
98  allocate(this%r(n))
99  allocate(this%ml(n))
100  allocate(this%mu(n))
101  allocate(this%wk1(n))
102 
103  allocate(this%c(this%lgmres))
104  allocate(this%s(this%lgmres))
105  allocate(this%gam(this%lgmres + 1))
106 
107  allocate(this%z(n,this%lgmres))
108  allocate(this%v(n,this%lgmres))
109 
110  allocate(this%h(this%lgmres,this%lgmres))
111 
112 
113  if (present(rel_tol) .and. present(abs_tol) .and. present(monitor)) then
114  call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
115  else if (present(rel_tol) .and. present(abs_tol)) then
116  call this%ksp_init(max_iter, rel_tol, abs_tol)
117  else if (present(monitor) .and. present(abs_tol)) then
118  call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
119  else if (present(rel_tol) .and. present(monitor)) then
120  call this%ksp_init(max_iter, rel_tol, monitor = monitor)
121  else if (present(rel_tol)) then
122  call this%ksp_init(max_iter, rel_tol = rel_tol)
123  else if (present(abs_tol)) then
124  call this%ksp_init(max_iter, abs_tol = abs_tol)
125  else if (present(monitor)) then
126  call this%ksp_init(max_iter, monitor = monitor)
127  else
128  call this%ksp_init(max_iter)
129  end if
130 
131  end subroutine sx_gmres_init
132 
134  subroutine sx_gmres_free(this)
135  class(sx_gmres_t), intent(inout) :: this
136 
137  call this%ksp_free()
138 
139  if (allocated(this%w)) then
140  deallocate(this%w)
141  end if
142 
143  if (allocated(this%c)) then
144  deallocate(this%c)
145  end if
146 
147  if (allocated(this%r)) then
148  deallocate(this%r)
149  end if
150 
151  if (allocated(this%z)) then
152  deallocate(this%z)
153  end if
154 
155  if (allocated(this%h)) then
156  deallocate(this%h)
157  end if
158 
159  if (allocated(this%ml)) then
160  deallocate(this%ml)
161  end if
162 
163  if (allocated(this%v)) then
164  deallocate(this%v)
165  end if
166 
167  if (allocated(this%s)) then
168  deallocate(this%s)
169  end if
170 
171  if (allocated(this%mu)) then
172  deallocate(this%mu)
173  end if
174 
175  if (allocated(this%gam)) then
176  deallocate(this%gam)
177  end if
178 
179  if (allocated(this%wk1)) then
180  deallocate(this%wk1)
181  end if
182 
183  nullify(this%M)
184 
185  end subroutine sx_gmres_free
186 
188  function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_results)
189  class(sx_gmres_t), intent(inout) :: this
190  class(ax_t), intent(in) :: ax
191  type(field_t), intent(inout) :: x
192  integer, intent(in) :: n
193  real(kind=rp), dimension(n), intent(in) :: f
194  type(coef_t), intent(inout) :: coef
195  type(bc_list_t), intent(in) :: blst
196  type(gs_t), intent(inout) :: gs_h
197  type(ksp_monitor_t) :: ksp_results
198  integer, optional, intent(in) :: niter
199  integer :: iter, max_iter, glb_n
200  integer :: i, j, k, ierr
201  real(kind=rp), parameter :: one = 1.0
202  real(kind=rp) :: rnorm
203  real(kind=rp) :: alpha, temp, l
204  real(kind=rp) :: ratio, div0, norm_fac
205  logical :: conv
206  integer outer
207 
208  conv = .false.
209  iter = 0
210  glb_n = n / x%msh%nelv * x%msh%glb_nelv
211 
212  if (present(niter)) then
213  max_iter = niter
214  else
215  max_iter = this%max_iter
216  end if
217 
218  call rone(this%ml, n)
219  call rone(this%mu, n)
220  norm_fac = one / sqrt(coef%volume)
221  call rzero(x%x, n)
222  call rzero(this%gam, this%lgmres + 1)
223  call rone(this%s, this%lgmres)
224  call rone(this%c, this%lgmres)
225  call rzero(this%h, this%lgmres * this%lgmres)
226  outer = 0
227  call this%monitor_start('GMRES')
228  do while (.not. conv .and. iter .lt. max_iter)
229  outer = outer + 1
230 
231  if(iter.eq.0) then
232  call col3(this%r,this%ml,f,n)
233  else
234  !update residual
235  call copy (this%r,f,n)
236  call ax%compute(this%w, x%x, coef, x%msh, x%Xh)
237  call gs_h%op(this%w, n, gs_op_add)
238  call bc_list_apply(blst, this%w, n)
239  call add2s2(this%r,this%w,-one,n)
240  call col2(this%r,this%ml,n)
241  endif
242  this%gam(1) = sqrt(glsc3(this%r, this%r, coef%mult, n))
243  if(iter.eq.0) then
244  div0 = this%gam(1) * norm_fac
245  ksp_results%res_start = div0
246  endif
247 
248  if (abscmp(this%gam(1), 0.0_rp)) return
249 
250  rnorm = 0.0_rp
251  temp = one / this%gam(1)
252  call cmult2(this%v(1,1), this%r, temp, n)
253  do j = 1, this%lgmres
254  iter = iter+1
255  call col3(this%w, this%mu, this%v(1,j), n)
256 
257  !Apply precond
258  call this%M%solve(this%z(1,j), this%w, n)
259 
260  call ax%compute(this%w, this%z(1,j), coef, x%msh, x%Xh)
261  call gs_h%op(this%w, n, gs_op_add)
262  call bc_list_apply(blst, this%w, n)
263  call col2(this%w, this%ml, n)
264 
265  do i = 1, j
266  this%h(i,j) = 0.0_rp
267  do k = 1, n
268  this%h(i,j) = this%h(i,j) + &
269  this%w(k) * this%v(k,i) * coef%mult(k,1,1,1)
270  end do
271  end do
272 
273  !Could probably be done inplace...
274  call mpi_allreduce(this%h(1,j), this%wk1, j, &
275  mpi_real_precision, mpi_sum, neko_comm, ierr)
276  call copy(this%h(1,j), this%wk1, j)
277 
278  do i = 1, j
279  do k = 1, n
280  this%w(k) = this%w(k) - this%h(i,j) * this%v(k,i)
281  end do
282  end do
283 
284  !apply Givens rotations to new column
285  do i=1,j-1
286  temp = this%h(i,j)
287  this%h(i ,j) = this%c(i)*temp + this%s(i)*this%h(i+1,j)
288  this%h(i+1,j) = -this%s(i)*temp + this%c(i)*this%h(i+1,j)
289  end do
290 
291  alpha = sqrt(glsc3(this%w, this%w, coef%mult, n))
292  rnorm = 0.0_rp
293  if(abscmp(alpha, 0.0_rp)) then
294  conv = .true.
295  exit
296  end if
297  l = sqrt(this%h(j,j) * this%h(j,j) + alpha**2)
298  temp = one / l
299  this%c(j) = this%h(j,j) * temp
300  this%s(j) = alpha * temp
301  this%h(j,j) = l
302  this%gam(j+1) = -this%s(j) * this%gam(j)
303  this%gam(j) = this%c(j) * this%gam(j)
304 
305  rnorm = abs(this%gam(j+1)) * norm_fac
306  call this%monitor_iter(iter, rnorm)
307  ratio = rnorm / div0
308  if (rnorm .lt. this%abs_tol) then
309  conv = .true.
310  exit
311  end if
312 
313  if (iter + 1 .gt. max_iter) exit
314 
315  if( j .lt. this%lgmres) then
316  temp = one / alpha
317  call cmult2(this%v(1,j+1), this%w, temp, n)
318  endif
319  end do
320  j = min(j, this%lgmres)
321  !back substitution
322  do k = j, 1, -1
323  temp = this%gam(k)
324  do i = j, k+1, -1
325  temp = temp - this%h(k,i) * this%c(i)
326  enddo
327  this%c(k) = temp / this%h(k,k)
328  enddo
329  !sum up Arnoldi vectors
330  do i = 1, j
331  do k = 1, n
332  x%x(k,1,1,1) = x%x(k,1,1,1) + this%c(i) * this%z(k,i)
333  end do
334  end do
335  end do
336  call this%monitor_stop()
337  ksp_results%res_final = rnorm
338  ksp_results%iter = iter
339  end function sx_gmres_solve
340 
342  function sx_gmres_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
343  n, coef, blstx, blsty, blstz, gs_h, niter) result(ksp_results)
344  class(sx_gmres_t), intent(inout) :: this
345  class(ax_t), intent(in) :: ax
346  type(field_t), intent(inout) :: x
347  type(field_t), intent(inout) :: y
348  type(field_t), intent(inout) :: z
349  integer, intent(in) :: n
350  real(kind=rp), dimension(n), intent(in) :: fx
351  real(kind=rp), dimension(n), intent(in) :: fy
352  real(kind=rp), dimension(n), intent(in) :: fz
353  type(coef_t), intent(inout) :: coef
354  type(bc_list_t), intent(in) :: blstx
355  type(bc_list_t), intent(in) :: blsty
356  type(bc_list_t), intent(in) :: blstz
357  type(gs_t), intent(inout) :: gs_h
358  type(ksp_monitor_t), dimension(3) :: ksp_results
359  integer, optional, intent(in) :: niter
360 
361  ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
362  ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
363  ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
364 
365  end function sx_gmres_solve_coupled
366 
367 end module gmres_sx
368 
369 
Defines a Matrix-vector product.
Definition: ax.f90:34
Defines a boundary condition.
Definition: bc.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_real_precision
MPI type for working precision of REAL types.
Definition: comm.F90:23
Defines a field.
Definition: field.f90:34
Gather-scatter.
Defines various GMRES methods.
Definition: gmres_sx.f90:34
subroutine sx_gmres_free(this)
Deallocate a standard GMRES solver.
Definition: gmres_sx.f90:135
subroutine sx_gmres_init(this, n, max_iter, M, lgmres, rel_tol, abs_tol, monitor)
Initialise a standard GMRES solver.
Definition: gmres_sx.f90:75
type(ksp_monitor_t) function, dimension(3) sx_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_sx.f90:344
type(ksp_monitor_t) function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
Standard PCG solve.
Definition: gmres_sx.f90:189
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:701
real(kind=rp) function, public glsc3(a, b, c, n)
Weighted inner product .
Definition: math.f90:895
subroutine, public rone(a, n)
Set all elements to one.
Definition: math.f90:228
subroutine, public col2(a, b, n)
Vector multiplication .
Definition: math.f90:729
subroutine, public copy(a, b, n)
Copy a vector .
Definition: math.f90:239
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition: math.f90:742
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
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 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
Standard preconditioned generalized minimal residual method (SX version)
Definition: gmres_sx.f90:49
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