66 use precon,
only :
pc_t, precon_factory, precon_destroy
86 krylov_solver_factory, krylov_solver_destroy
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)
then
156 call neko_error(
'Polynomial order < 2 not supported for hsmg precon')
163 allocate(this%grids(this%nlvls))
164 allocate(this%w(dof%size()))
165 allocate(this%r(dof%size()))
169 call msh%all_deformed()
172 call this%e%init(dof,
'work array')
173 call this%wf%init(dof,
'work 2')
175 call this%Xh_crs%init(
gll, lx_crs, lx_crs, lx_crs)
176 call this%dm_crs%init(msh, this%Xh_crs)
177 call this%gs_crs%init(this%dm_crs)
178 call this%e_crs%init(this%dm_crs,
'work crs')
179 call this%c_crs%init(this%gs_crs)
181 call this%Xh_mg%init(
gll, lx_mid, lx_mid, lx_mid)
182 call this%dm_mg%init(msh, this%Xh_mg)
183 call this%gs_mg%init(this%dm_mg)
184 call this%e_mg%init(this%dm_mg,
'work midl')
185 call this%c_mg%init(this%gs_mg)
188 call ax_helm_factory(this%ax, full_formulation = .false.)
191 call precon_factory(this%pc_crs,
'jacobi')
194 if (
present(crs_pctype))
then
195 call krylov_solver_factory(this%crs_solver, &
196 this%dm_crs%size(), trim(crs_pctype),
ksp_max_iter, m = this%pc_crs)
198 call krylov_solver_factory(this%crs_solver, &
199 this%dm_crs%size(),
'cg',
ksp_max_iter, m = this%pc_crs)
202 call this%bc_crs%init_base(this%c_crs)
203 call this%bc_mg%init_base(this%c_mg)
204 call this%bc_reg%init_base(coef)
205 if (bclst%n .gt. 0)
then
207 call this%bc_reg%mark_facets(bclst%bc(i)%bcp%marked_facet)
208 call this%bc_crs%mark_facets(bclst%bc(i)%bcp%marked_facet)
209 call this%bc_mg%mark_facets(bclst%bc(i)%bcp%marked_facet)
212 call this%bc_reg%finalize()
213 call this%bc_reg%set_g(
real(0d0,
rp))
217 call this%bc_crs%finalize()
218 call this%bc_crs%set_g(
real(0d0,
rp))
223 call this%bc_mg%finalize()
224 call this%bc_mg%set_g(0.0_rp)
228 call this%schwarz%init(xh, dof, gs_h, this%bclst_reg, msh)
229 call this%schwarz_mg%init(this%Xh_mg, this%dm_mg, this%gs_mg,&
232 call this%interp_fine_mid%init(xh, this%Xh_mg)
233 call this%interp_mid_crs%init(this%Xh_mg, this%Xh_crs)
235 call hsmg_fill_grid(dof, gs_h, xh, coef, this%bclst_reg, this%schwarz, &
236 this%e, this%grids, 3)
237 call hsmg_fill_grid(this%dm_mg, this%gs_mg, this%Xh_mg, this%c_mg, &
238 this%bclst_mg, this%schwarz_mg, this%e_mg, &
241 this%c_crs, this%bclst_crs, this%schwarz_crs, &
242 this%e_crs, this%grids, 1)
250 select type(pc => this%pc_crs)
252 call pc%init(this%c_crs, this%dm_crs, this%gs_crs)
254 call pc%init(this%c_crs, this%dm_crs, this%gs_crs)
256 call pc%init(this%c_crs, this%dm_crs, this%gs_crs)
264 class(
hsmg_t),
intent(inout) :: this
267 this%grids(1)%coef%ifh2 = .false.
268 call copy(this%grids(1)%coef%h1, this%grids(3)%coef%h1, &
269 this%grids(1)%dof%size())
271 call device_copy(this%grids(1)%coef%h1_d, this%grids(3)%coef%h1_d, &
272 this%grids(1)%dof%size())
278 type(
dofmap_t),
target,
intent(in) :: dof
279 type(
gs_t),
target,
intent(in) :: gs_h
280 type(
space_t),
target,
intent(in) :: Xh
281 type(
coef_t),
target,
intent(in) :: coef
282 type(
bc_list_t),
target,
intent(in) :: bclst
283 type(
schwarz_t),
target,
intent(in) :: schwarz
284 type(
field_t),
target,
intent(in) :: e
285 integer,
intent(in) :: l
286 type(
multigrid_t),
intent(inout),
dimension(l) :: grids
290 grids(l)%gs_h => gs_h
292 grids(l)%coef => coef
293 grids(l)%bclst => bclst
300 class(
hsmg_t),
intent(inout) :: this
302 if (
allocated(this%ax))
then
306 if (
allocated(this%grids))
then
307 deallocate(this%grids)
310 if (
allocated(this%w))
then
314 if (
allocated(this%r))
then
318 call this%schwarz%free()
319 call this%schwarz_mg%free()
321 call this%c_crs%free()
322 call this%c_mg%free()
324 call this%e_mg%free()
325 call this%e_crs%free()
327 call this%gs_crs%free()
328 call this%gs_mg%free()
329 call this%interp_mid_crs%free()
330 call this%interp_fine_mid%free()
332 if (
allocated(this%crs_solver))
then
333 call krylov_solver_destroy(this%crs_solver)
334 deallocate(this%crs_solver)
337 if (
allocated(this%pc_crs))
then
338 call precon_destroy(this%pc_crs)
345 integer,
intent(in) :: n
346 class(
hsmg_t),
intent(inout) :: this
347 real(kind=
rp),
dimension(n),
intent(inout) :: z
348 real(kind=
rp),
dimension(n),
intent(inout) :: r
349 type(c_ptr) :: z_d, r_d
351 integer :: thrdid, nthrds
363 call device_col2(this%r_d, this%grids(3)%coef%mult_d, &
364 this%grids(3)%dof%size())
366 call this%interp_fine_mid%map(this%e%x, this%r, &
367 this%msh%nelv, this%grids(2)%Xh)
368 call this%grids(2)%gs_h%op(this%e%x, &
369 this%grids(2)%dof%size(), gs_op_add, this%gs_event)
373 call device_copy(this%w_d, this%e%x_d, this%grids(2)%dof%size())
375 this%grids(2)%dof%size())
377 call device_col2(this%w_d, this%grids(2)%coef%mult_d, &
378 this%grids(2)%dof%size())
380 call this%interp_mid_crs%map(this%wf%x, this%w, this%msh%nelv, &
383 call device_copy(this%w_d, this%e%x_d, this%grids(2)%dof%size())
385 this%grids(2)%dof%size())
394 if (thrdid .eq. 0)
then
396 call this%grids(3)%schwarz%compute(z, this%r)
397 call this%grids(2)%schwarz%compute(this%grids(2)%e%x, this%w)
400 if (nthrds .eq. 1 .or. thrdid .eq. 1)
then
402 call this%grids(1)%gs_h%op(this%wf%x, &
403 this%grids(1)%dof%size(), gs_op_add, this%gs_event)
406 this%grids(1)%dof%size())
408 crs_info = this%crs_solver%solve(this%Ax, this%grids(1)%e, &
410 this%grids(1)%dof%size(), &
411 this%grids(1)%coef, &
412 this%grids(1)%bclst, &
413 this%grids(1)%gs_h, this%niter)
416 this%grids(1)%dof%size())
421 call this%interp_mid_crs%map(this%w, this%grids(1)%e%x, &
422 this%msh%nelv, this%grids(2)%Xh)
423 call device_add2(this%grids(2)%e%x_d, this%w_d, this%grids(2)%dof%size())
425 call this%interp_fine_mid%map(this%w, this%grids(2)%e%x, &
426 this%msh%nelv, this%grids(3)%Xh)
427 call device_add2(z_d, this%w_d, this%grids(3)%dof%size())
428 call this%grids(3)%gs_h%op(z, this%grids(3)%dof%size(), &
429 gs_op_add, this%gs_event)
432 this%grids(3)%dof%size())
435 call copy(this%r, r, n)
438 call this%grids(3)%schwarz%compute(z, this%r)
440 call col2(this%r, this%grids(3)%coef%mult, &
441 this%grids(3)%dof%size())
443 call this%interp_fine_mid%map(this%w, this%r, &
444 this%msh%nelv, this%grids(2)%Xh)
445 call this%grids(2)%gs_h%op(this%w, this%grids(2)%dof%size(), gs_op_add)
447 call this%grids(2)%schwarz%compute(this%grids(2)%e%x, this%w)
448 call col2(this%w, this%grids(2)%coef%mult, this%grids(2)%dof%size())
450 call this%interp_mid_crs%map(this%r, this%w, &
451 this%msh%nelv, this%grids(1)%Xh)
454 call this%grids(1)%gs_h%op(this%r, this%grids(1)%dof%size(), gs_op_add)
456 this%grids(1)%dof%size())
458 crs_info = this%crs_solver%solve(this%Ax, this%grids(1)%e, this%r, &
459 this%grids(1)%dof%size(), &
460 this%grids(1)%coef, &
461 this%grids(1)%bclst, &
462 this%grids(1)%gs_h, this%niter)
465 this%grids(1)%dof%size())
468 call this%interp_mid_crs%map(this%w, this%grids(1)%e%x, &
469 this%msh%nelv, this%grids(2)%Xh)
470 call add2(this%grids(2)%e%x, this%w, this%grids(2)%dof%size())
472 call this%interp_fine_mid%map(this%w, this%grids(2)%e%x, &
473 this%msh%nelv, this%grids(3)%Xh)
474 call add2(z, this%w, this%grids(3)%dof%size())
475 call this%grids(3)%gs_h%op(z, this%grids(3)%dof%size(), gs_op_add)
476 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)
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)
Vector addition .
subroutine, public device_col2(a_d, b_d, n)
Vector multiplication .
subroutine, public device_copy(a_d, b_d, n)
Copy a vector .
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.
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_device
integer, parameter, public rp
Global precision used in computations.
subroutine, public profiler_start_region(name, region_id)
Started a named (name) profiler region.
subroutine, public profiler_end_region(name, region_id)
End the most recently started 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.
Base type for a boundary condition.
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.