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 logical :: zero_initial_guess = .false.
71 subroutine cheby_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
72 class(
cheby_t),
intent(inout),
target :: this
73 integer,
intent(in) :: max_iter
74 class(
pc_t),
optional,
intent(in),
target :: M
75 integer,
intent(in) :: n
76 real(kind=
rp),
optional,
intent(in) :: rel_tol
77 real(kind=
rp),
optional,
intent(in) :: abs_tol
78 logical,
optional,
intent(in) :: monitor
89 if (
present(rel_tol) .and.
present(abs_tol) .and.
present(monitor))
then
90 call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
91 else if (
present(rel_tol) .and.
present(abs_tol))
then
92 call this%ksp_init(max_iter, rel_tol, abs_tol)
93 else if (
present(monitor) .and.
present(abs_tol))
then
94 call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
95 else if (
present(rel_tol) .and.
present(monitor))
then
96 call this%ksp_init(max_iter, rel_tol, monitor = monitor)
97 else if (
present(rel_tol))
then
98 call this%ksp_init(max_iter, rel_tol = rel_tol)
99 else if (
present(abs_tol))
then
100 call this%ksp_init(max_iter, abs_tol = abs_tol)
101 else if (
present(monitor))
then
102 call this%ksp_init(max_iter, monitor = monitor)
104 call this%ksp_init(max_iter)
110 class(
cheby_t),
intent(inout) :: this
111 if (
allocated(this%d))
then
117 class(
cheby_t),
intent(inout) :: this
118 class(
ax_t),
intent(in) :: Ax
119 type(
field_t),
intent(inout) :: x
120 integer,
intent(in) :: n
121 type(
coef_t),
intent(inout) :: coef
123 type(
gs_t),
intent(inout) :: gs_h
124 real(kind=
rp) :: lam, b, a, rn
125 real(kind=
rp) :: boost = 1.1_rp
126 real(kind=
rp) :: lam_factor = 30.0_rp
127 real(kind=
rp) :: wtw, dtw, dtd
129 associate(w => this%w, d => this%d, r => this%r)
133 call random_number(rn)
136 call gs_h%op(d, n, gs_op_add)
137 call blst%apply(d, n)
140 do i = 1, this%power_its
141 call ax%compute(w, d, coef, x%msh, x%Xh)
142 call gs_h%op(w, n, gs_op_add)
143 call blst%apply(w, n)
144 if (
associated(this%schwarz))
then
145 call this%schwarz%compute(r, w)
148 call this%M%solve(r, w, n)
152 wtw =
glsc3(w, coef%mult, w, n)
153 call cmult2(d, w, 1.0_rp/sqrt(wtw), n)
154 call blst%apply(d, n)
157 call ax%compute(w, d, coef, x%msh, x%Xh)
158 call gs_h%op(w, n, gs_op_add)
159 call blst%apply(w, n)
160 if (
associated(this%schwarz))
then
161 call this%schwarz%compute(r, w)
164 call this%M%solve(r, w, n)
168 dtw =
glsc3(d, coef%mult, w, n)
169 dtd =
glsc3(d, coef%mult, d, n)
173 this%tha = (b+a)/2.0_rp
174 this%dlt = (b-a)/2.0_rp
176 this%recompute_eigs = .false.
181 function cheby_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
183 class(
cheby_t),
intent(inout) :: this
184 class(ax_t),
intent(in) :: ax
185 type(field_t),
intent(inout) :: x
186 integer,
intent(in) :: n
187 real(kind=rp),
dimension(n),
intent(in) :: f
188 type(coef_t),
intent(inout) :: coef
189 type(bc_list_t),
intent(inout) :: blst
190 type(gs_t),
intent(inout) :: gs_h
191 type(ksp_monitor_t) :: ksp_results
192 integer,
optional,
intent(in) :: niter
193 integer :: iter, max_iter
194 real(kind=rp) :: a, b, rtr, rnorm, norm_fac
196 if (this%recompute_eigs)
then
200 if (
present(niter))
then
203 max_iter = this%max_iter
205 norm_fac = 1.0_rp / sqrt(coef%volume)
207 associate( w => this%w, r => this%r, d => this%d)
210 call ax%compute(w, x%x, coef, x%msh, x%Xh)
211 call gs_h%op(w, n, gs_op_add)
212 call blst%apply(w, n)
215 rtr = glsc3(r, coef%mult, r, n)
216 rnorm = sqrt(rtr) * norm_fac
217 ksp_results%res_start = rnorm
218 ksp_results%res_final = rnorm
222 call this%M%solve(w, r, n)
224 a = 2.0_rp / this%tha
225 call add2s2(x%x, d, a, n)
228 do iter = 2, max_iter
231 call ax%compute(w, x%x, coef, x%msh, x%Xh)
232 call gs_h%op(w, n, gs_op_add)
233 call blst%apply(w, n)
236 call this%M%solve(w, r, n)
238 if (iter .eq. 2)
then
239 b = 0.5_rp * (this%dlt * a)**2
241 b = (this%dlt * a / 2.0_rp)**2
243 a = 1.0_rp/(this%tha - b/a)
244 call add2s1(d, w, b, n)
246 call add2s2(x%x, d, a, n)
251 call ax%compute(w, x%x, coef, x%msh, x%Xh)
252 call gs_h%op(w, n, gs_op_add)
253 call blst%apply(w, n)
255 rtr = glsc3(r, coef%mult, r, n)
256 rnorm = sqrt(rtr) * norm_fac
257 ksp_results%res_final = rnorm
258 ksp_results%iter = iter
259 ksp_results%converged = this%is_converged(iter, rnorm)
264 function cheby_impl(this, Ax, x, f, n, coef, blst, gs_h, niter) &
266 class(
cheby_t),
intent(inout) :: this
267 class(ax_t),
intent(in) :: ax
268 type(field_t),
intent(inout) :: x
269 integer,
intent(in) :: n
270 real(kind=rp),
dimension(n),
intent(in) :: f
271 type(coef_t),
intent(inout) :: coef
272 type(bc_list_t),
intent(inout) :: blst
273 type(gs_t),
intent(inout) :: gs_h
274 type(ksp_monitor_t) :: ksp_results
275 integer,
optional,
intent(in) :: niter
276 integer :: iter, max_iter
277 real(kind=rp) :: a, b, rtr, rnorm, norm_fac
278 real(kind=rp) :: rhok, rhokp1, sig1, tmp1, tmp2
280 if (this%recompute_eigs)
then
284 if (
present(niter))
then
287 max_iter = this%max_iter
289 norm_fac = 1.0_rp / sqrt(coef%volume)
291 associate( w => this%w, r => this%r, d => this%d)
293 if (.not.this%zero_initial_guess)
then
294 call ax%compute(w, x%x, coef, x%msh, x%Xh)
295 call gs_h%op(w, n, gs_op_add)
296 call blst%apply(w, n)
297 call sub3(r, f, w, n)
300 this%zero_initial_guess = .false.
304 if (
associated(this%schwarz))
then
305 call this%schwarz%compute(d, r)
307 call this%M%solve(d, r, n)
309 call cmult( d, (1.0_rp / this%tha), n)
310 call add2( x%x, d, n)
312 sig1 = this%tha / this%dlt
316 do iter = 2, max_iter
317 rhokp1 = 1.0_rp / (2.0_rp * sig1 - rhok)
319 tmp2 = 2.0_rp * rhokp1 / this%dlt
322 call ax%compute(w, x%x, coef, x%msh, x%Xh)
323 call gs_h%op(w, n, gs_op_add)
324 call blst%apply(w, n)
325 call sub3(r, f, w, n)
327 if (
associated(this%schwarz))
then
328 call this%schwarz%compute(w, r)
330 call this%M%solve(w, r, n)
332 call cmult( d, tmp1, n)
333 call add2s2( d, w, tmp2, n)
334 call add2( x%x, d, n)
342 n, coef, blstx, blsty, blstz, gs_h, niter)
result(ksp_results)
343 class(
cheby_t),
intent(inout) :: this
344 class(ax_t),
intent(in) :: ax
345 type(field_t),
intent(inout) :: x
346 type(field_t),
intent(inout) :: y
347 type(field_t),
intent(inout) :: z
348 integer,
intent(in) :: n
349 real(kind=rp),
dimension(n),
intent(in) :: fx
350 real(kind=rp),
dimension(n),
intent(in) :: fy
351 real(kind=rp),
dimension(n),
intent(in) :: fz
352 type(coef_t),
intent(inout) :: coef
353 type(bc_list_t),
intent(inout) :: blstx
354 type(bc_list_t),
intent(inout) :: blsty
355 type(bc_list_t),
intent(inout) :: blstz
356 type(gs_t),
intent(inout) :: gs_h
357 type(ksp_monitor_t),
dimension(3) :: ksp_results
358 integer,
optional,
intent(in) :: niter
360 ksp_results(1) = this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
361 ksp_results(2) = this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
362 ksp_results(3) = this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
__device__ T solve(const T u, const T y, const T guess, const T nu, const T kappa, const T B)
Defines a Matrix-vector product.
Chebyshev preconditioner.
type(ksp_monitor_t) function cheby_impl(this, ax, x, f, n, coef, blst, gs_h, niter)
A chebyshev preconditioner.
subroutine cheby_free(this)
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.
subroutine cheby_init(this, n, max_iter, m, rel_tol, abs_tol, monitor)
Initialise a standard solver.
subroutine cheby_power(this, ax, x, n, coef, blst, gs_h)
type(ksp_monitor_t) function cheby_solve(this, ax, x, f, n, coef, blst, gs_h, niter)
A chebyshev preconditioner.
Implements the base abstract type for Krylov solvers plus helper types.
subroutine, public cmult(a, c, n)
Multiplication by constant c .
subroutine, public cmult2(a, b, c, n)
Multiplication by constant c .
real(kind=rp) function, public glsc3(a, b, c, n)
Weighted inner product .
subroutine, public add2s1(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
real(kind=rp) function, public glsc2(a, b, n)
Weighted inner product .
subroutine, public rone(a, n)
Set all elements to one.
subroutine, public sub3(a, b, c, n)
Vector subtraction .
subroutine, public add2(a, b, n)
Vector addition .
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
subroutine, public sub2(a, b, n)
Vector substraction .
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
integer, parameter, public rp
Global precision used in computations.
Overlapping schwarz solves.
Defines a function space.
Base type for a matrix-vector product providing .
A list of allocatable `bc_t`. Follows the standard interface of lists.
Defines a Chebyshev preconditioner.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Type for storing initial and final residuals in a Krylov solver.
Base abstract type for a canonical Krylov method, solving .
Defines a canonical Krylov preconditioner.
The function space for the SEM solution fields.