Neko  0.9.0
A portable framework for high-order spectral element flow simulations
cheby.f90
Go to the documentation of this file.
1 ! Copyright (c) 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 cheby
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 mesh, only : mesh_t
42  use space, only : space_t
43  use gather_scatter, only : gs_t, gs_op_add
44  use bc, only : bc_list_t, bc_list_apply
45  use math, only : glsc3, rzero, rone, copy, sub2, cmult2, abscmp, glsc2, &
46  add2s1, add2s2
47  use comm
48  implicit none
49  private
50 
52  type, public, extends(ksp_t) :: cheby_t
53  real(kind=rp), allocatable :: d(:)
54  real(kind=rp), allocatable :: w(:)
55  real(kind=rp), allocatable :: r(:)
56  real(kind=rp) :: tha, dlt
57  integer :: power_its = 150
58  logical :: recompute_eigs = .true.
59  contains
60  procedure, pass(this) :: init => cheby_init
61  procedure, pass(this) :: free => cheby_free
62  procedure, pass(this) :: solve => cheby_solve
63  procedure, pass(this) :: solve_coupled => cheby_solve_coupled
64  end type cheby_t
65 
66 contains
67 
69  subroutine cheby_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
70  class(cheby_t), intent(inout), target :: this
71  integer, intent(in) :: max_iter
72  class(pc_t), optional, intent(inout), target :: M
73  integer, intent(in) :: n
74  real(kind=rp), optional, intent(inout) :: rel_tol
75  real(kind=rp), optional, intent(inout) :: abs_tol
76  logical, optional, intent(in) :: monitor
77 
78  call this%free()
79  allocate(this%d(n))
80  allocate(this%w(n))
81  allocate(this%r(n))
82 
83  if (present(m)) then
84  this%M => m
85  end if
86 
87  if (present(rel_tol) .and. present(abs_tol) .and. present(monitor)) then
88  call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
89  else if (present(rel_tol) .and. present(abs_tol)) then
90  call this%ksp_init(max_iter, rel_tol, abs_tol)
91  else if (present(monitor) .and. present(abs_tol)) then
92  call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
93  else if (present(rel_tol) .and. present(monitor)) then
94  call this%ksp_init(max_iter, rel_tol, monitor = monitor)
95  else if (present(rel_tol)) then
96  call this%ksp_init(max_iter, rel_tol = rel_tol)
97  else if (present(abs_tol)) then
98  call this%ksp_init(max_iter, abs_tol = abs_tol)
99  else if (present(monitor)) then
100  call this%ksp_init(max_iter, monitor = monitor)
101  else
102  call this%ksp_init(max_iter)
103  end if
104 
105  end subroutine cheby_init
106 
107  subroutine cheby_free(this)
108  class(cheby_t), intent(inout) :: this
109  if (allocated(this%d)) then
110  deallocate(this%d)
111  end if
112  end subroutine cheby_free
113 
114  subroutine cheby_power(this, Ax, x, n, coef, blst, gs_h)
115  class(cheby_t), intent(inout) :: this
116  class(ax_t), intent(inout) :: Ax
117  type(field_t), intent(inout) :: x
118  integer, intent(in) :: n
119  type(coef_t), intent(inout) :: coef
120  type(bc_list_t), intent(inout) :: blst
121  type(gs_t), intent(inout) :: gs_h
122  real(kind=rp) :: lam, b, a, rn
123  real(kind=rp) :: boost = 1.2_rp
124  real(kind=rp) :: lam_factor = 30.0_rp
125  real(kind=rp) :: wtw, dtw, dtd
126  integer :: i
127  associate(w => this%w, d => this%d)
128 
129  do i = 1, n
130  !TODO: replace with a better way to initialize power method
131  call random_number(rn)
132  d(i) = rn + 10.0_rp
133  end do
134  call gs_h%op(d, n, gs_op_add)
135  call bc_list_apply(blst, d, n)
136 
137  !Power method to get lamba max
138  do i = 1, this%power_its
139  call ax%compute(w, d, coef, x%msh, x%Xh)
140  call gs_h%op(w, n, gs_op_add)
141  call bc_list_apply(blst, w, n)
142 
143  wtw = glsc3(w, coef%mult, w, n)
144  call cmult2(d, w, 1.0_rp/sqrt(wtw), n)
145  call bc_list_apply(blst, d, n)
146  end do
147 
148  call ax%compute(w, d, coef, x%msh, x%Xh)
149  call gs_h%op(w, n, gs_op_add)
150  call bc_list_apply(blst, w, n)
151 
152  dtw = glsc3(d, coef%mult, w, n)
153  dtd = glsc3(d, coef%mult, d, n)
154  lam = dtw / dtd
155  b = lam * boost
156  a = lam / lam_factor
157  this%tha = (b+a)/2.0_rp
158  this%dlt = (b-a)/2.0_rp
159 
160  this%recompute_eigs = .false.
161  end associate
162  end subroutine cheby_power
163 
165  function cheby_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
166  result(ksp_results)
167  class(cheby_t), intent(inout) :: this
168  class(ax_t), intent(inout) :: ax
169  type(field_t), intent(inout) :: x
170  integer, intent(in) :: n
171  real(kind=rp), dimension(n), intent(inout) :: f
172  type(coef_t), intent(inout) :: coef
173  type(bc_list_t), intent(inout) :: blst
174  type(gs_t), intent(inout) :: gs_h
175  type(ksp_monitor_t) :: ksp_results
176  integer, optional, intent(in) :: niter
177  integer :: iter, max_iter
178  real(kind=rp) :: a, b, rtr, rnorm, norm_fac
179 
180  if (this%recompute_eigs) then
181  call cheby_power(this, ax, x, n, coef, blst, gs_h)
182  end if
183 
184  if (present(niter)) then
185  max_iter = niter
186  else
187  max_iter = this%max_iter
188  end if
189  norm_fac = 1.0_rp / sqrt(coef%volume)
190 
191  associate( w => this%w, r => this%r, d => this%d)
192  ! calculate residual
193  call copy(r, f, n)
194  call ax%compute(w, x%x, coef, x%msh, x%Xh)
195  call gs_h%op(w, n, gs_op_add)
196  call bc_list_apply(blst, w, n)
197  call sub2(r, w, n)
198 
199  rtr = glsc3(r, coef%mult, r, n)
200  rnorm = sqrt(rtr) * norm_fac
201  ksp_results%res_start = rnorm
202  ksp_results%res_final = rnorm
203  ksp_results%iter = 0
204 
205  ! First iteration
206  call this%M%solve(w, r, n)
207  call copy(d, w, n)
208  a = 2.0_rp / this%tha
209  call add2s2(x%x, d, a, n)! x = x + a*d
210 
211  ! Rest of the iterations
212  do iter = 2, max_iter
213  ! calculate residual
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 bc_list_apply(blst, w, n)
218  call sub2(r, w, n)
219 
220  call this%M%solve(w, r, n)
221 
222  if (iter .eq. 2) then
223  b = 0.5_rp * (this%dlt * a)**2
224  else
225  b = (this%dlt * a / 2.0_rp)**2
226  end if
227  a = 1.0_rp/(this%tha - b/a)
228  call add2s1(d, w, b, n)! d = w + b*d
229 
230  call add2s2(x%x, d, a, n)! x = x + a*d
231  end do
232 
233  ! calculate residual
234  call copy(r, f, n)
235  call ax%compute(w, x%x, coef, x%msh, x%Xh)
236  call gs_h%op(w, n, gs_op_add)
237  call bc_list_apply(blst, w, n)
238  call sub2(r, w, n)
239  rtr = glsc3(r, coef%mult, r, n)
240  rnorm = sqrt(rtr) * norm_fac
241  ksp_results%res_final = rnorm
242  ksp_results%iter = iter
243  end associate
244  end function cheby_solve
245 
247  function cheby_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
248  n, coef, blstx, blsty, blstz, gs_h, niter) result(ksp_results)
249  class(cheby_t), intent(inout) :: this
250  class(ax_t), intent(inout) :: ax
251  type(field_t), intent(inout) :: x
252  type(field_t), intent(inout) :: y
253  type(field_t), intent(inout) :: z
254  integer, intent(in) :: n
255  real(kind=rp), dimension(n), intent(inout) :: fx
256  real(kind=rp), dimension(n), intent(inout) :: fy
257  real(kind=rp), dimension(n), intent(inout) :: fz
258  type(coef_t), intent(inout) :: coef
259  type(bc_list_t), intent(inout) :: blstx
260  type(bc_list_t), intent(inout) :: blsty
261  type(bc_list_t), intent(inout) :: blstz
262  type(gs_t), intent(inout) :: gs_h
263  type(ksp_monitor_t), dimension(3) :: ksp_results
264  integer, optional, intent(in) :: niter
265 
266  ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
267  ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
268  ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
269 
270  end function cheby_solve_coupled
271 
272 end module cheby
273 
274 
Defines a Matrix-vector product.
Definition: ax.f90:34
Defines a boundary condition.
Definition: bc.f90:34
Chebyshev preconditioner.
Definition: cheby.f90:34
subroutine cheby_power(this, Ax, x, n, coef, blst, gs_h)
Definition: cheby.f90:115
subroutine cheby_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
Initialise a standard solver.
Definition: cheby.f90:70
type(ksp_monitor_t) function cheby_solve(this, Ax, x, f, n, coef, blst, gs_h, niter)
A chebyshev preconditioner.
Definition: cheby.f90:167
subroutine cheby_free(this)
Definition: cheby.f90:108
type(ksp_monitor_t) function, dimension(3) cheby_solve_coupled(this, Ax, x, y, z, fx, fy, fz, n, coef, blstx, blsty, blstz, gs_h, niter)
Standard Chebyshev coupled solve.
Definition: cheby.f90:249
Coefficients.
Definition: coef.f90:34
Definition: comm.F90:1
Defines a field.
Definition: field.f90:34
Gather-scatter.
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 add2s1(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
Definition: math.f90:658
real(kind=rp) function, public glsc2(a, b, n)
Weighted inner product .
Definition: math.f90:876
subroutine, public rone(a, n)
Set all elements to one.
Definition: math.f90:228
subroutine, public copy(a, b, n)
Copy a vector .
Definition: math.f90:239
subroutine, public rzero(a, n)
Zero a real vector.
Definition: math.f90:195
subroutine, public sub2(a, b, n)
Vector substraction .
Definition: math.f90:629
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
Definition: math.f90:673
Defines a mesh.
Definition: mesh.f90:34
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Krylov preconditioner.
Definition: precon.f90:34
Defines a function space.
Definition: space.f90:34
Base type for a matrix-vector product providing .
Definition: ax.f90:43
A list of boundary conditions.
Definition: bc.f90:104
Defines a Chebyshev preconditioner.
Definition: cheby.f90:52
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
The function space for the SEM solution fields.
Definition: space.f90:62