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.
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
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)
102 call this%ksp_init(max_iter)
108 class(
cheby_t),
intent(inout) :: this
109 if (
allocated(this%d))
then
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
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
127 associate(w => this%w, d => this%d)
131 call random_number(rn)
134 call gs_h%op(d, n, gs_op_add)
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)
143 wtw =
glsc3(w, coef%mult, w, n)
144 call cmult2(d, w, 1.0_rp/sqrt(wtw), n)
148 call ax%compute(w, d, coef, x%msh, x%Xh)
149 call gs_h%op(w, n, gs_op_add)
152 dtw =
glsc3(d, coef%mult, w, n)
153 dtd =
glsc3(d, coef%mult, d, n)
157 this%tha = (b+a)/2.0_rp
158 this%dlt = (b-a)/2.0_rp
160 this%recompute_eigs = .false.
165 function cheby_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
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
180 if (this%recompute_eigs)
then
184 if (
present(niter))
then
187 max_iter = this%max_iter
189 norm_fac = 1.0_rp / sqrt(coef%volume)
191 associate( w => this%w, r => this%r, d => this%d)
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)
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
206 call this%M%solve(w, r, n)
208 a = 2.0_rp / this%tha
209 call add2s2(x%x, d, a, n)
212 do iter = 2, max_iter
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)
220 call this%M%solve(w, r, n)
222 if (iter .eq. 2)
then
223 b = 0.5_rp * (this%dlt * a)**2
225 b = (this%dlt * a / 2.0_rp)**2
227 a = 1.0_rp/(this%tha - b/a)
228 call add2s1(d, w, b, n)
230 call add2s2(x%x, d, a, 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)
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
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
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)
Defines a Matrix-vector product.
Defines a boundary condition.
Chebyshev preconditioner.
subroutine cheby_power(this, Ax, x, n, coef, blst, gs_h)
subroutine cheby_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
Initialise a standard solver.
type(ksp_monitor_t) function cheby_solve(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.
Implements the base abstract type for Krylov solvers plus helper types.
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 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.
Defines a function space.
Base type for a matrix-vector product providing .
A list of boundary conditions.
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.