Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
tree_amg_multigrid.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!
45 use num_types, only: rp
46 use utils, only : neko_error
47 use math, only : add2, rzero, glsc2, sub3, col2
50 use comm
51 use coefs, only : coef_t
52 use mesh, only : mesh_t
53 use space, only : space_t
54 use ax_product, only: ax_t
55 use bc_list, only : bc_list_t
56 use gather_scatter, only : gs_t, gs_op_add
61 use logger, only : neko_log, log_size
65 use, intrinsic :: iso_c_binding
66 implicit none
67 private
68
69 type :: tamg_wrk_t
70 real(kind=rp), allocatable :: r(:)
71 real(kind=rp), allocatable :: rc(:)
72 real(kind=rp), allocatable :: tmp(:)
73 type(c_ptr) :: r_d = c_null_ptr
74 type(c_ptr) :: rc_d = c_null_ptr
75 type(c_ptr) :: tmp_d = c_null_ptr
76 end type tamg_wrk_t
77
79 type, public :: tamg_solver_t
80 type(tamg_hierarchy_t), allocatable :: amg
81 type(amg_cheby_t), allocatable :: smoo(:)
82 type(tamg_wrk_t), allocatable :: wrk(:)
83 !type(amg_jacobi_t), allocatable :: jsmoo(:)
84 integer :: nlvls
85 integer :: max_iter
86 contains
87 procedure, pass(this) :: init => tamg_mg_init
88 procedure, pass(this) :: solve => tamg_mg_solve
89 end type tamg_solver_t
90
91contains
92
102 subroutine tamg_mg_init(this, ax, Xh, coef, msh, gs_h, nlvls, blst, &
103 max_iter, cheby_degree)
104 class(tamg_solver_t), intent(inout), target :: this
105 class(ax_t), target, intent(in) :: ax
106 type(space_t), target, intent(in) :: Xh
107 type(coef_t), target, intent(in) :: coef
108 type(mesh_t), target, intent(in) :: msh
109 type(gs_t), target, intent(in) :: gs_h
110 type(bc_list_t), target, intent(in) :: blst
111 integer, intent(in) :: nlvls
112 integer, intent(in) :: max_iter
113 integer, intent(in) :: cheby_degree
114 integer :: lvl, n, mlvl, target_num_aggs
115 integer, allocatable :: agg_nhbr(:,:), asdf(:,:)
116 character(len=LOG_SIZE) :: log_buf
117
118 call neko_log%section('AMG')
119
120 write(log_buf, '(A28,I2,A8)') 'Creating AMG hierarchy with', &
121 nlvls, 'levels.'
122 call neko_log%message(log_buf)
123
124 allocate( this%amg )
125 call this%amg%init(ax, xh, coef, msh, gs_h, nlvls, blst)
126
128 call aggregate_finest_level(this%amg, xh%lx, xh%ly, xh%lz, msh%nelv)
129
131 allocate( agg_nhbr, source = msh%facet_neigh )
132 do mlvl = 2, nlvls-1
133 target_num_aggs = this%amg%lvl(mlvl-1)%nnodes / 8
134 call print_preagg_info( mlvl, target_num_aggs)
135 if ( target_num_aggs .lt. 4 ) then
136 call neko_error( &
137 "TAMG: Too many levels. Not enough DOFs for coarsest grid.")
138 end if
139 call aggregate_greedy( this%amg, mlvl, target_num_aggs, agg_nhbr, asdf)
140 agg_nhbr = asdf
141 deallocate( asdf )
142 end do
143 deallocate( agg_nhbr )
144
146 call aggregate_end(this%amg, nlvls)
147
148 this%max_iter = max_iter
149
150 this%nlvls = this%amg%nlvls!TODO: read from parameter
151 if (this%nlvls .gt. this%amg%nlvls) then
152 call neko_error( &
153 "Requested number multigrid levels &
154 & is greater than the initialized AMG levels")
155 end if
156
157 allocate(this%smoo(0:(nlvls)))
158 do lvl = 0, nlvls-1
159 n = this%amg%lvl(lvl+1)%fine_lvl_dofs
160 call this%smoo(lvl)%init(n, lvl, cheby_degree)
161 end do
162
164 allocate(this%wrk(nlvls))
165 do lvl = 1, nlvls
166 n = this%amg%lvl(lvl)%fine_lvl_dofs
167 allocate( this%wrk(lvl)%r(n) )
168 allocate( this%wrk(lvl)%rc(n) )
169 allocate( this%wrk(lvl)%tmp(n) )
170 if (neko_bcknd_device .eq. 1) then
171 call device_map( this%wrk(lvl)%r, this%wrk(lvl)%r_d, n)
172 call device_map( this%wrk(lvl)%rc, this%wrk(lvl)%rc_d, n)
173 call device_map( this%wrk(lvl)%tmp, this%wrk(lvl)%tmp_d, n)
174 end if
175 end do
176
177 !allocate(this%jsmoo(0:(nlvls)))
178 !do lvl = 0, nlvls-1
179 ! n = this%amg%lvl(lvl+1)%fine_lvl_dofs
180 ! call this%jsmoo(lvl)%init(n ,lvl, cheby_degree)
181 !end do
182
183 call fill_lvl_map(this%amg)
184
185 call neko_log%end_section()
186
187 end subroutine tamg_mg_init
188
189
194 subroutine tamg_mg_solve(this, z, r, n)
195 integer, intent(in) :: n
196 class(tamg_solver_t), intent(inout) :: this
197 real(kind=rp), dimension(n), intent(inout) :: z
198 real(kind=rp), dimension(n), intent(inout) :: r
199 type(c_ptr) :: z_d
200 type(c_ptr) :: r_d
201 integer :: iter, max_iter
202 logical :: zero_initial_guess
203
204 max_iter = this%max_iter
205
206 if (neko_bcknd_device .eq. 1) then
207 z_d = device_get_ptr(z)
208 r_d = device_get_ptr(r)
209 ! Zero out the initial guess becuase we do not handle null spaces very well...
210 call device_rzero(z_d, n)
211 zero_initial_guess = .true.
212 ! Call the amg cycle
213 do iter = 1, max_iter
214 call tamg_mg_cycle_d(z, r, z_d, r_d, n, 0, this%amg, this, &
215 zero_initial_guess)
216 zero_initial_guess = .false.
217 end do
218 else
219 ! Zero out the initial guess becuase we do not handle null spaces very well...
220 call rzero(z, n)
221 zero_initial_guess = .true.
222 ! Call the amg cycle
223 do iter = 1, max_iter
224 call tamg_mg_cycle(z, r, n, 0, this%amg, this, &
225 zero_initial_guess)
226 zero_initial_guess = .false.
227 end do
228 end if
229 end subroutine tamg_mg_solve
230
231
239 recursive subroutine tamg_mg_cycle(x, b, n, lvl, amg, mgstuff, &
240 zero_initial_guess)
241 integer, intent(in) :: n
242 real(kind=rp), intent(inout) :: x(n)
243 real(kind=rp), intent(inout) :: b(n)
244 type(tamg_hierarchy_t), intent(inout) :: amg
245 type(tamg_solver_t), intent(inout) :: mgstuff
246 logical, intent(in) :: zero_initial_guess
247 integer, intent(in) :: lvl
248 real(kind=rp) :: r(n)
249 real(kind=rp) :: rc(n)
250 real(kind=rp) :: tmp(n)
251 integer :: iter, num_iter
252 integer :: max_lvl
253 integer :: i, cyt
254 max_lvl = mgstuff%nlvls-1
258 call mgstuff%smoo(lvl)%solve(x, b, n, amg, &
259 zero_initial_guess)
260 if (lvl .eq. max_lvl) then
261 return
262 end if
266 call calc_resid(r, x, b, amg, lvl, n)
270 call amg%interp_f2c(rc, r, lvl+1)
274 call rzero(tmp, n)
275 call tamg_mg_cycle(tmp, rc, amg%lvl(lvl+1)%nnodes, lvl+1, amg, mgstuff, &
276 .true.)
280 call amg%interp_c2f(r, tmp, lvl+1)
284 call add2(x, r, n)
288 call mgstuff%smoo(lvl)%solve(x,b, n, amg)
289 end subroutine tamg_mg_cycle
290
298 recursive subroutine tamg_mg_cycle_d(x, b, x_d, b_d, n, lvl, amg, mgstuff, &
299 zero_initial_guess)
300 integer, intent(in) :: n
301 real(kind=rp), intent(inout) :: x(n)
302 real(kind=rp), intent(inout) :: b(n)
303 type(c_ptr) :: x_d
304 type(c_ptr) :: b_d
305 type(tamg_hierarchy_t), intent(inout) :: amg
306 type(tamg_solver_t), intent(inout) :: mgstuff
307 logical, intent(in) :: zero_initial_guess
308 integer, intent(in) :: lvl
309 integer :: iter, num_iter
310 integer :: max_lvl
311 integer :: i, cyt
312 max_lvl = mgstuff%nlvls-1
316 call mgstuff%smoo(lvl)%device_solve(x, b, x_d, b_d, n, amg, &
317 zero_initial_guess)
318 if (lvl .eq. max_lvl) then
319 return
320 end if
321 associate( r => mgstuff%wrk(lvl+1)%r, r_d => mgstuff%wrk(lvl+1)%r_d, &
322 rc => mgstuff%wrk(lvl+1)%rc, rc_d => mgstuff%wrk(lvl+1)%rc_d, &
323 tmp => mgstuff%wrk(lvl+1)%tmp, tmp_d => mgstuff%wrk(lvl+1)%tmp_d )
327 call amg%device_matvec(r, x, r_d, x_d, lvl)
328 call device_sub3(r_d, b_d, r_d, n)
332 call amg%interp_f2c_d(rc_d, r_d, lvl+1)
336 call device_rzero(tmp_d, n)
337 call tamg_mg_cycle_d(tmp, rc, tmp_d, rc_d, &
338 amg%lvl(lvl+1)%nnodes, lvl+1, amg, mgstuff, .true.)
342 call amg%interp_c2f_d(r_d, tmp_d, lvl+1, r)
346 call device_add2(x_d, r_d, n)
350 call mgstuff%smoo(lvl)%device_solve(x, b, x_d, b_d, n, amg)
351 end associate
352 end subroutine tamg_mg_cycle_d
353
354
362 subroutine calc_resid(r, x, b, amg, lvl, n)
363 integer, intent(in) :: n
364 real(kind=rp), intent(inout) :: r(n)
365 real(kind=rp), intent(inout) :: x(n)
366 real(kind=rp), intent(inout) :: b(n)
367 type(tamg_hierarchy_t), intent(inout) :: amg
368 integer, intent(in) :: lvl
369 integer :: i
370 call amg%matvec(r, x, lvl)
371 call sub3(r, b, r, n)
372 end subroutine calc_resid
373
374
375 subroutine print_preagg_info(lvl, nagg)
376 integer, intent(in) :: lvl, nagg
377 character(len=LOG_SIZE) :: log_buf
378 !TODO: calculate min and max agg size
379 write(log_buf, '(A8,I2,A31)') '-- level', lvl, &
380 '-- Calling Greedy Aggregation'
381 call neko_log%message(log_buf)
382 write(log_buf, '(A33,I6)') 'Target Aggregates:', nagg
383 call neko_log%message(log_buf)
384 end subroutine print_preagg_info
385
386 subroutine print_resid_info(r, x, b, r_d, x_d, b_d, amg, lvl, n)
387 integer, intent(in) :: lvl, n
388 real(kind=rp), intent(inout) :: r(n)
389 real(kind=rp), intent(inout) :: x(n)
390 real(kind=rp), intent(inout) :: b(n)
391 type(c_ptr) :: r_d
392 type(c_ptr) :: x_d
393 type(c_ptr) :: b_d
394 type(tamg_hierarchy_t), intent(inout) :: amg
395 real(kind=rp) :: val
396 character(len=LOG_SIZE) :: log_buf
397
398 call amg%device_matvec(r, x, r_d, x_d, lvl)
399 call device_sub3(r_d, b_d, r_d, n)
400 val = device_glsc2(r_d, r_d, n)
401
402 write(log_buf, '(A33,I6,F12.6)') 'tAMG resid:', lvl, val
403 call neko_log%message(log_buf)
404 end subroutine print_resid_info
405
406 subroutine fill_lvl_map(amg)
407 type(tamg_hierarchy_t), intent(inout) :: amg
408 integer :: i, j, k, l, nid, n
409 do j = 1, amg%lvl(1)%nnodes
410 do k = 1, amg%lvl(1)%nodes(j)%ndofs
411 nid = amg%lvl(1)%nodes(j)%dofs(k)
412 amg%lvl(1)%map_finest2lvl(nid) = amg%lvl(1)%nodes(j)%gid
413 end do
414 end do
415 n = size(amg%lvl(1)%map_finest2lvl)
416 do l = 2, amg%nlvls
417 do i = 1, n
418 nid = amg%lvl(l-1)%map_finest2lvl(i)
419 do j = 1, amg%lvl(l)%nnodes
420 do k = 1, amg%lvl(l)%nodes(j)%ndofs
421 if (nid .eq. amg%lvl(l)%nodes(j)%dofs(k)) then
422 amg%lvl(l)%map_finest2lvl(i) = amg%lvl(l)%nodes(j)%gid
423 end if
424 end do
425 end do
426 end do
427 end do
428 if (neko_bcknd_device .eq. 1) then
429 do l = 1, amg%nlvls
430 amg%lvl(l)%map_finest2lvl(0) = n
431 call device_memcpy( amg%lvl(l)%map_finest2lvl, &
432 amg%lvl(l)%map_finest2lvl_d, n, &
433 host_to_device, .true.)
434 call device_memcpy( amg%lvl(l)%map_f2c, &
435 amg%lvl(l)%map_f2c_d, amg%lvl(l)%fine_lvl_dofs+1, &
436 host_to_device, .true.)
437 end do
438 end if
439 end subroutine fill_lvl_map
440end module tree_amg_multigrid
__device__ T solve(const T u, const T y, const T guess, const T nu, const T kappa, const T B)
Return the device pointer for an associated Fortran array.
Definition device.F90:96
Map a Fortran array to a device (allocate and associate)
Definition device.F90:72
Copy data between host and device (or device and device)
Definition device.F90:66
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
subroutine, public device_sub3(a_d, b_d, c_d, n, strm)
Vector subtraction .
subroutine, public device_add2(a_d, b_d, n, strm)
Vector addition .
subroutine, public device_rzero(a_d, n, strm)
Zero a real vector.
subroutine, public device_col2(a_d, b_d, n, strm)
Vector multiplication .
real(kind=rp) function, public device_glsc2(a_d, b_d, n, strm)
Weighted inner product .
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:47
subroutine, public device_free(x_d)
Deallocate memory on the device.
Definition device.F90:214
Gather-scatter.
Logging routines.
Definition log.f90:34
type(log_t), public neko_log
Global log stream.
Definition log.f90:70
integer, parameter, public log_size
Definition log.f90:46
Definition math.f90:60
real(kind=rp) function, public glsc2(a, b, n)
Weighted inner product .
Definition math.f90:1007
subroutine, public sub3(a, b, c, n)
Vector subtraction .
Definition math.f90:787
subroutine, public add2(a, b, n)
Vector addition .
Definition math.f90:732
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:860
subroutine, public rzero(a, n)
Zero a real vector.
Definition math.f90:211
Defines a mesh.
Definition mesh.f90:34
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Defines a function space.
Definition space.f90:34
Implements an aggregation for TreeAMG hierarchy structure.
subroutine aggregate_end(tamg, lvl_id)
Aggregate all dofs to a single point to form a tree-like structure.
subroutine aggregate_greedy(tamg, lvl_id, max_aggs, facet_neigh, agg_nhbr)
Aggregates dofs based on adjacent dofs.
subroutine aggregate_finest_level(tamg, lx, ly, lz, ne)
Aggregaiton on finest level Aggregates all dofs in an element into a single aggregate.
Implements multigrid using the TreeAMG hierarchy structure. USE:
recursive subroutine tamg_mg_cycle(x, b, n, lvl, amg, mgstuff, zero_initial_guess)
Recrsive multigrid cycle for the TreeAMG solver object.
subroutine calc_resid(r, x, b, amg, lvl, n)
Wrapper function to calculate residyal.
subroutine print_preagg_info(lvl, nagg)
recursive subroutine tamg_mg_cycle_d(x, b, x_d, b_d, n, lvl, amg, mgstuff, zero_initial_guess)
Recrsive multigrid cycle for the TreeAMG solver object on device.
subroutine fill_lvl_map(amg)
subroutine print_resid_info(r, x, b, r_d, x_d, b_d, amg, lvl, n)
subroutine tamg_mg_solve(this, z, r, n)
Solver function for the TreeAMG solver object.
subroutine tamg_mg_init(this, ax, xh, coef, msh, gs_h, nlvls, blst, max_iter, cheby_degree)
Initialization of the TreeAMG multigrid solver.
Implements smoothers for use with TreeAMG matrix vector product.
Implements the base type for TreeAMG hierarchy structure.
Definition tree_amg.f90:34
subroutine, public tamg_node_init(node, gid, ndofs)
Initialization of a TreeAMG tree node.
Definition tree_amg.f90:203
subroutine, public tamg_lvl_init(tamg_lvl, lvl, nnodes, ndofs)
Initialization of a TreeAMG level.
Definition tree_amg.f90:171
Utilities.
Definition utils.f90:35
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:48
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
The function space for the SEM solution fields.
Definition space.f90:62
Type for a TreeAMG hierarchy.
Definition tree_amg.f90:94
Type for the TreeAMG solver.
Type for Chebyshev iteration using TreeAMG matvec.