94 type(
gs_t),
pointer :: gs_h
115 class(
ksp_t),
allocatable :: crs_solver
116 integer :: niter = 10
117 class(
pc_t),
allocatable :: pc_crs
119 real(kind=
rp),
allocatable :: r(:)
122 real(kind=
rp),
allocatable :: w(:)
123 type(c_ptr) :: w_d = c_null_ptr
124 type(c_ptr) :: r_d = c_null_ptr
125 type(c_ptr) :: hsmg_event
126 type(c_ptr) :: gs_event
137 subroutine hsmg_init(this, msh, Xh, coef, dof, gs_h, bclst, crs_pctype)
138 class(
hsmg_t),
intent(inout),
target :: this
139 type(
mesh_t),
intent(inout),
target :: msh
140 type(
space_t),
intent(inout),
target :: Xh
141 type(
coef_t),
intent(inout),
target :: coef
142 type(
dofmap_t),
intent(inout),
target :: dof
143 type(
gs_t),
intent(inout),
target :: gs_h
144 type(
bc_list_t),
intent(inout),
target :: bclst
145 character(len=*),
optional :: crs_pctype
147 integer :: lx_crs, lx_mid
152 if (xh%lx .lt. 5)
then
153 lx_mid =
max(xh%lx-1,3)
155 if(xh%lx .le. 2)
call neko_error(
'Polynomial order < 2 not supported for hsmg precon')
161 allocate(this%grids(this%nlvls))
162 allocate(this%w(dof%size()))
163 allocate(this%r(dof%size()))
167 call msh%all_deformed()
170 call this%e%init(dof,
'work array')
171 call this%wf%init(dof,
'work 2')
173 call this%Xh_crs%init(
gll, lx_crs, lx_crs, lx_crs)
174 this%dm_crs =
dofmap_t(msh, this%Xh_crs)
175 call this%gs_crs%init(this%dm_crs)
176 call this%e_crs%init(this%dm_crs,
'work crs')
177 call this%c_crs%init(this%gs_crs)
179 call this%Xh_mg%init(
gll, lx_mid, lx_mid, lx_mid)
180 this%dm_mg =
dofmap_t(msh, this%Xh_mg)
181 call this%gs_mg%init(this%dm_mg)
182 call this%e_mg%init(this%dm_mg,
'work midl')
183 call this%c_mg%init(this%gs_mg)
202 if (
present(crs_pctype))
then
204 this%dm_crs%size(), trim(crs_pctype),
ksp_max_iter, m = this%pc_crs)
207 this%dm_crs%size(),
'cg',
ksp_max_iter, m = this%pc_crs)
210 call this%bc_crs%init(this%c_crs)
211 call this%bc_mg%init(this%c_mg)
212 call this%bc_reg%init(coef)
213 if (bclst%n .gt. 0)
then
215 call this%bc_reg%mark_facets(bclst%bc(i)%bcp%marked_facet)
216 call this%bc_crs%mark_facets(bclst%bc(i)%bcp%marked_facet)
217 call this%bc_mg%mark_facets(bclst%bc(i)%bcp%marked_facet)
220 call this%bc_reg%finalize()
221 call this%bc_reg%set_g(
real(0d0,
rp))
225 call this%bc_crs%finalize()
226 call this%bc_crs%set_g(
real(0d0,
rp))
231 call this%bc_mg%finalize()
232 call this%bc_mg%set_g(0.0_rp)
236 call this%schwarz%init(xh, dof, gs_h, this%bclst_reg, msh)
237 call this%schwarz_mg%init(this%Xh_mg, this%dm_mg, this%gs_mg,&
240 call this%interp_fine_mid%init(xh,this%Xh_mg)
241 call this%interp_mid_crs%init(this%Xh_mg,this%Xh_crs)
243 call hsmg_fill_grid(dof, gs_h, xh, coef, this%bclst_reg, this%schwarz, &
244 this%e, this%grids, 3)
245 call hsmg_fill_grid(this%dm_mg, this%gs_mg, this%Xh_mg, this%c_mg, &
246 this%bclst_mg, this%schwarz_mg, this%e_mg, &
249 this%c_crs, this%bclst_crs, this%schwarz_crs, &
250 this%e_crs, this%grids, 1)
257 select type(pc => this%pc_crs)
259 call pc%init(this%c_crs, this%dm_crs, this%gs_crs)
261 call pc%init(this%c_crs, this%dm_crs, this%gs_crs)
263 call pc%init(this%c_crs, this%dm_crs, this%gs_crs)
270 class(
hsmg_t),
intent(inout) :: this
273 this%grids(1)%coef%ifh2 = .false.
274 call copy(this%grids(1)%coef%h1, this%grids(3)%coef%h1, &
275 this%grids(1)%dof%size())
277 call device_copy(this%grids(1)%coef%h1_d, this%grids(3)%coef%h1_d, &
278 this%grids(1)%dof%size())
284 type(
dofmap_t),
target,
intent(in):: dof
285 type(
gs_t),
target,
intent(in) :: gs_h
286 type(
space_t),
target,
intent(in) :: Xh
287 type(
coef_t),
target,
intent(in) :: coef
288 type(
bc_list_t),
target,
intent(in) :: bclst
289 type(
schwarz_t),
target,
intent(in) :: schwarz
290 type(
field_t),
target,
intent(in) :: e
291 integer,
intent(in) :: l
292 type(
multigrid_t),
intent(inout),
dimension(l) :: grids
296 grids(l)%gs_h => gs_h
298 grids(l)%coef => coef
299 grids(l)%bclst => bclst
306 class(
hsmg_t),
intent(inout) :: this
308 if (
allocated(this%ax))
then
312 if (
allocated(this%grids))
then
313 deallocate(this%grids)
316 if (
allocated(this%w))
then
320 if (
allocated(this%r))
then
324 call this%schwarz%free()
325 call this%schwarz_mg%free()
327 call this%c_crs%free()
328 call this%c_mg%free()
330 call this%e_mg%free()
331 call this%e_crs%free()
333 call this%gs_crs%free()
334 call this%gs_mg%free()
335 call this%interp_mid_crs%free()
336 call this%interp_fine_mid%free()
338 if (
allocated(this%crs_solver))
then
340 deallocate(this%crs_solver)
343 if (
allocated(this%pc_crs))
then
344 select type(pc => this%pc_crs)
350 deallocate(this%pc_crs)
357 integer,
intent(in) :: n
358 class(
hsmg_t),
intent(inout) :: this
359 real(kind=
rp),
dimension(n),
intent(inout) :: z
360 real(kind=
rp),
dimension(n),
intent(inout) :: r
361 type(c_ptr) :: z_d, r_d
363 integer :: i, thrdid, nthrds
375 call device_col2(this%r_d, this%grids(3)%coef%mult_d, &
376 this%grids(3)%dof%size())
378 call this%interp_fine_mid%map(this%e%x, this%r, &
379 this%msh%nelv, this%grids(2)%Xh)
380 call this%grids(2)%gs_h%op(this%e%x, &
381 this%grids(2)%dof%size(), gs_op_add, this%gs_event)
385 call device_copy(this%w_d, this%e%x_d, this%grids(2)%dof%size())
387 this%grids(2)%dof%size())
389 call device_col2(this%w_d, this%grids(2)%coef%mult_d, &
390 this%grids(2)%dof%size())
392 call this%interp_mid_crs%map(this%wf%x, this%w,this%msh%nelv, &
395 call device_copy(this%w_d, this%e%x_d, this%grids(2)%dof%size())
397 this%grids(2)%dof%size())
406 if (thrdid .eq. 0)
then
408 call this%grids(3)%schwarz%compute(z, this%r)
409 call this%grids(2)%schwarz%compute(this%grids(2)%e%x,this%w)
412 if (nthrds .eq. 1 .or. thrdid .eq. 1)
then
414 call this%grids(1)%gs_h%op(this%wf%x, &
415 this%grids(1)%dof%size(), gs_op_add, this%gs_event)
418 this%grids(1)%dof%size())
420 crs_info = this%crs_solver%solve(this%Ax, this%grids(1)%e, this%wf%x, &
421 this%grids(1)%dof%size(), &
422 this%grids(1)%coef, &
423 this%grids(1)%bclst, &
424 this%grids(1)%gs_h, this%niter)
427 this%grids(1)%dof%size())
432 call this%interp_mid_crs%map(this%w, this%grids(1)%e%x, &
433 this%msh%nelv, this%grids(2)%Xh)
434 call device_add2(this%grids(2)%e%x_d, this%w_d, this%grids(2)%dof%size())
436 call this%interp_fine_mid%map(this%w, this%grids(2)%e%x, &
437 this%msh%nelv, this%grids(3)%Xh)
438 call device_add2(z_d, this%w_d, this%grids(3)%dof%size())
439 call this%grids(3)%gs_h%op(z, this%grids(3)%dof%size(), &
440 gs_op_add, this%gs_event)
442 call device_col2(z_d, this%grids(3)%coef%mult_d, this%grids(3)%dof%size())
445 call copy(this%r, r, n)
448 call this%grids(3)%schwarz%compute(z, this%r)
450 call col2(this%r, this%grids(3)%coef%mult, &
451 this%grids(3)%dof%size())
453 call this%interp_fine_mid%map(this%w, this%r, &
454 this%msh%nelv, this%grids(2)%Xh)
455 call this%grids(2)%gs_h%op(this%w, this%grids(2)%dof%size(), gs_op_add)
457 call this%grids(2)%schwarz%compute(this%grids(2)%e%x,this%w)
458 call col2(this%w, this%grids(2)%coef%mult, this%grids(2)%dof%size())
460 call this%interp_mid_crs%map(this%r,this%w,this%msh%nelv,this%grids(1)%Xh)
463 call this%grids(1)%gs_h%op(this%r, this%grids(1)%dof%size(), gs_op_add)
465 this%grids(1)%dof%size())
467 crs_info = this%crs_solver%solve(this%Ax, this%grids(1)%e, this%r, &
468 this%grids(1)%dof%size(), &
469 this%grids(1)%coef, &
470 this%grids(1)%bclst, &
471 this%grids(1)%gs_h, this%niter)
474 this%grids(1)%dof%size())
477 call this%interp_mid_crs%map(this%w, this%grids(1)%e%x, &
478 this%msh%nelv, this%grids(2)%Xh)
479 call add2(this%grids(2)%e%x, this%w, this%grids(2)%dof%size())
481 call this%interp_fine_mid%map(this%w, this%grids(2)%e%x, &
482 this%msh%nelv, this%grids(3)%Xh)
483 call add2(z, this%w, this%grids(3)%dof%size())
484 call this%grids(3)%gs_h%op(z, this%grids(3)%dof%size(), gs_op_add)
485 call col2(z, this%grids(3)%coef%mult, this%grids(3)%dof%size())
Return the device pointer for an associated Fortran array.
Map a Fortran array to a device (allocate and associate)
subroutine, public ax_helm_factory(Ax)
Factory routine for the a Helmholtz problem matrix-vector product. The selection is based on the comp...
Defines a Matrix-vector product.
Defines a boundary condition.
subroutine, public bc_list_add(bclst, bc)
Add a condition to a list of boundary conditions.
subroutine, public bc_list_init(bclst, size)
Constructor for a list of boundary conditions.
subroutine, public bc_list_apply_scalar(bclst, x, n, t, tstep)
Apply a list of boundary conditions to a scalar field.
Jacobi preconditioner accelerator backend.
subroutine, public device_add2(a_d, b_d, n)
subroutine, public device_col2(a_d, b_d, n)
subroutine, public device_copy(a_d, b_d, n)
Device abstraction, common interface for various accelerators.
subroutine, public device_event_sync(event)
Synchronize an event.
subroutine, public device_event_create(event, flags)
Create a device event queue.
Defines a dirichlet boundary condition.
Defines a mapping of the degrees of freedom.
subroutine hsmg_solve(this, z, r, n)
The h1mg preconditioner from Nek5000.
subroutine hsmg_init(this, msh, Xh, coef, dof, gs_h, bclst, crs_pctype)
subroutine hsmg_set_h(this)
subroutine hsmg_free(this)
subroutine hsmg_fill_grid(dof, gs_h, Xh, coef, bclst, schwarz, e, grids, l)
Routines to interpolate between different spaces.
subroutine, public krylov_solver_destroy(ksp)
Destroy an interative Krylov solver.
subroutine, public krylov_solver_factory(ksp, n, solver, max_iter, abstol, M)
Initialize an iterative Krylov solver.
Implements the base abstract type for Krylov solvers plus helper types.
integer, parameter, public ksp_max_iter
Maximum number of iters.
subroutine, public add2(a, b, n)
Vector addition .
subroutine, public col2(a, b, n)
Vector multiplication .
subroutine, public copy(a, b, n)
Copy a vector .
integer, parameter neko_bcknd_sx
integer, parameter neko_bcknd_device
integer, parameter neko_bcknd_xsmm
integer, parameter, public rp
Global precision used in computations.
subroutine, public profiler_end_region
End the most recently started profiler region.
subroutine, public profiler_start_region(name, region_id)
Started a named (name) profiler region.
Overlapping schwarz solves.
Defines a function space.
integer, parameter, public gll
Jacobi preconditioner SX-Aurora backend.
Base type for a matrix-vector product providing .
A list of boundary conditions.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Defines a jacobi preconditioner.
Generic Dirichlet boundary condition on .
Interpolation between two space::space_t.
Defines a jacobi preconditioner.
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.
Defines a jacobi preconditioner for SX-Aurora.