77 integer,
allocatable :: imin(:)
78 integer,
allocatable :: imax(:)
79 integer,
allocatable :: ri(:, :)
85 type(stencil_t),
public :: stencil
86 type(mesh_t),
pointer :: mesh => null()
87 integer,
allocatable :: nn(:)
88 integer,
public :: np = 0
90 real(real64),
allocatable,
public :: w(:,:)
92 logical,
public :: const_w = .
true.
94 type(accel_mem_t),
public :: buff_weights
95 type(accel_mem_t),
public :: buff_half_weights
97 character(len=40) :: label
100 integer,
public :: nri = 0
101 integer,
allocatable,
public :: ri(:,:)
102 integer,
allocatable,
public :: rimap(:)
103 integer,
allocatable,
public :: rimap_inv(:)
105 integer :: ninner = 0
106 integer :: nouter = 0
108 type(nl_operator_index_t) :: inner
109 type(nl_operator_index_t) :: outer
111 type(accel_kernel_t) :: kernel
112 type(accel_mem_t) :: buff_imin
113 type(accel_mem_t) :: buff_imax
114 type(accel_mem_t) :: buff_ri
115 type(accel_mem_t) :: buff_map
116 type(accel_mem_t) :: buff_all
117 type(accel_mem_t) :: buff_inner
118 type(accel_mem_t) :: buff_outer
119 type(accel_mem_t) :: buff_stencil
120 type(accel_mem_t) :: buff_ip_to_xyz
121 type(accel_mem_t) :: buff_xyz_to_ip
124 type(nl_operator_t),
public,
pointer :: coarser => null()
128 integer,
parameter :: &
134 integer,
parameter :: &
144 integer,
intent(in) :: opid, type
148 integer :: dfunction_global = -1
149 integer :: zfunction_global = -1
150 integer :: function_opencl
158 type(namespace_t),
intent(in) :: namespace
236 character(len=*),
intent(in) :: label
264 safe_allocate_source_a(opo%nn, opi%nn)
265 safe_allocate_source_a(opo%w, opi%w)
267 opo%const_w = opi%const_w
270 assert(
allocated(opi%ri))
272 safe_allocate_source_a(opo%ri, opi%ri)
273 safe_allocate_source_a(opo%rimap, opi%rimap)
274 safe_allocate_source_a(opo%rimap_inv, opi%rimap_inv)
276 if (opi%mesh%parallel_in_domains)
then
277 opo%inner%nri = opi%inner%nri
278 safe_allocate_source_a(opo%inner%imin, opi%inner%imin)
279 safe_allocate_source_a(opo%inner%imax, opi%inner%imax)
280 safe_allocate_source_a(opo%inner%ri, opi%inner%ri)
282 opo%outer%nri = opi%outer%nri
283 safe_allocate_source_a(opo%outer%imin, opi%outer%imin)
284 safe_allocate_source_a(opo%outer%imax, opi%outer%imax)
285 safe_allocate_source_a(opo%outer%ri, opi%outer%ri)
303 class(
space_t),
intent(in) :: space
304 type(
mesh_t),
target,
intent(in) :: mesh
306 integer,
intent(in) :: np
307 logical,
optional,
intent(in) :: const_w
308 logical,
optional,
intent(in) :: regenerate
310 integer :: ii, jj, p1(space%dim), time, current
311 integer,
allocatable :: st1(:), st2(:), st1r(:)
313 integer :: ir, maxp, iinner, iouter
314 logical :: change, force_change
315 character(len=200) :: flags
316 integer,
allocatable :: inner_points(:), outer_points(:), all_points(:)
320 if (mesh%parallel_in_domains .and. .not. const_w)
then
330 if (
present(const_w)) op%const_w = const_w
334 safe_allocate(op%w(1:op%stencil%size, 1))
335 message(1) =
'Debug: nl_operator_build: working with constant weights.'
338 safe_allocate(op%w(1:op%stencil%size, 1:op%np))
339 message(1) =
'Debug: nl_operator_build: working with non-constant weights.'
347 safe_allocate(st1(1:op%stencil%size))
348 safe_allocate(st1r(1:op%stencil%size))
349 safe_allocate(st2(1:op%stencil%size))
358 do jj = 1, op%stencil%size
365 st1(1:op%stencil%size) = st1(1:op%stencil%size) - ii
367 change = any(st1 /= st2)
371 force_change = any(st1 + ii > mesh%np) .and. all(st2 + ii - 1 <= mesh%np)
374 if (change .or. force_change)
then
379 if (time == 1) op%nri = op%nri + 1
383 current = current + 1
384 op%ri(1:op%stencil%size, current) = st1(1:op%stencil%size)
388 if (time == 2) op%rimap(ii) = current
394 safe_deallocate_a(op%ri)
395 safe_deallocate_a(op%rimap)
396 safe_deallocate_a(op%rimap_inv)
398 safe_allocate(op%ri(1:op%stencil%size, 1:op%nri))
399 safe_allocate(op%rimap(1:op%np))
400 safe_allocate(op%rimap_inv(1:op%nri + 1))
407 if (mesh%use_curvilinear)
then
408 safe_allocate(op%nn(1:op%nri))
410 op%nn = op%stencil%size
419 op%rimap_inv(op%rimap(jj) + 1) = jj
421 op%rimap_inv(op%nri + 1) = op%np
423 safe_deallocate_a(st1)
424 safe_deallocate_a(st1r)
425 safe_deallocate_a(st2)
428 nn = op%rimap_inv(jj + 1) - op%rimap_inv(jj)
431 if (op%mesh%parallel_in_domains)
then
438 maxp = op%rimap_inv(ir + 1) + maxval(op%ri(1:op%stencil%size, ir))
441 op%inner%nri = op%inner%nri + 1
442 assert(op%inner%nri <= op%nri)
445 op%outer%nri = op%outer%nri + 1
446 assert(op%outer%nri <= op%nri)
450 assert(op%inner%nri + op%outer%nri == op%nri)
453 safe_deallocate_a(op%inner%imin)
454 safe_deallocate_a(op%inner%imax)
455 safe_deallocate_a(op%inner%ri)
456 safe_deallocate_a(op%outer%imin)
457 safe_deallocate_a(op%outer%imax)
458 safe_deallocate_a(op%outer%ri)
460 safe_allocate(op%inner%imin(1:op%inner%nri + 1))
461 safe_allocate(op%inner%imax(1:op%inner%nri))
462 safe_allocate(op%inner%ri(1:op%stencil%size, 1:op%inner%nri))
464 safe_allocate(op%outer%imin(1:op%outer%nri + 1))
465 safe_allocate(op%outer%imax(1:op%outer%nri))
466 safe_allocate(op%outer%ri(1:op%stencil%size, 1:op%outer%nri))
472 maxp = op%rimap_inv(ir + 1) + maxval(op%ri(1:op%stencil%size, ir))
476 op%inner%imin(iinner) = op%rimap_inv(ir)
477 op%inner%imax(iinner) = op%rimap_inv(ir + 1)
478 op%inner%ri(1:op%stencil%size, iinner) = op%ri(1:op%stencil%size, ir)
482 op%outer%imin(iouter) = op%rimap_inv(ir)
483 op%outer%imax(iouter) = op%rimap_inv(ir + 1)
484 op%outer%ri(1:op%stencil%size, iouter) = op%ri(1:op%stencil%size, ir)
489 do ir = 1, op%inner%nri
490 do ii = op%inner%imin(ir) + 1, op%inner%imax(ir)
491 assert(all(ii + op%inner%ri(1:op%stencil%size, ir) <= mesh%np))
499 write(flags,
'(i5)') op%stencil%size
500 flags=
'-DNDIM=3 -DSTENCIL_SIZE='//trim(adjustl(flags))
502 if (op%mesh%parallel_in_domains) flags =
'-DINDIRECT '//trim(flags)
504 select case (function_opencl)
515 select case (function_opencl)
527 if (op%mesh%parallel_in_domains)
then
529 safe_allocate(inner_points(1:op%mesh%np))
530 safe_allocate(outer_points(1:op%mesh%np))
531 safe_allocate(all_points(1:op%mesh%np))
536 do ii = 1, op%mesh%np
537 all_points(ii) = ii - 1
538 maxp = ii + maxval(op%ri(1:op%stencil%size, op%rimap(ii)))
539 if (maxp <= op%mesh%np)
then
540 op%ninner = op%ninner + 1
541 inner_points(op%ninner) = ii - 1
543 op%nouter = op%nouter + 1
544 outer_points(op%nouter) = ii - 1
557 safe_deallocate_a(inner_points)
558 safe_deallocate_a(outer_points)
559 safe_deallocate_a(all_points)
573 integer :: istencil, idir
577 write(
message(1),
'(3a)')
'Debug info: Finite difference weights for ', trim(this%label),
'.'
578 write(
message(2),
'(a)')
' Spacing:'
579 do idir = 1, this%mesh%box%dim
580 write(
message(2),
'(a,f16.8)') trim(
message(2)), this%mesh%spacing(idir)
584 do istencil = 1, this%stencil%size
585 select case(this%mesh%box%dim)
587 write(
message(1),
'(a,i3,1i4,f25.10)')
' ', istencil, this%stencil%points(1:1, istencil), this%w(istencil, 1)
589 write(
message(1),
'(a,i3,2i4,f25.10)')
' ', istencil, this%stencil%points(1:2, istencil), this%w(istencil, 1)
591 write(
message(1),
'(a,i3,3i4,f25.10)')
' ', istencil, this%stencil%points(1:3, istencil), this%w(istencil, 1)
609 select case (function_opencl)
616 if (op%mesh%parallel_in_domains)
then
633 safe_deallocate_a(op%inner%imin)
634 safe_deallocate_a(op%inner%imax)
635 safe_deallocate_a(op%inner%ri)
636 safe_deallocate_a(op%outer%imin)
637 safe_deallocate_a(op%outer%imax)
638 safe_deallocate_a(op%outer%ri)
640 safe_deallocate_a(op%w)
642 safe_deallocate_a(op%ri)
643 safe_deallocate_a(op%rimap)
644 safe_deallocate_a(op%rimap_inv)
645 safe_deallocate_a(op%nn)
656 integer,
intent(in) :: is
657 integer,
intent(in) :: ip
659 res = ip + op%ri(is, op%rimap(ip))
670 if (accel_is_enabled() .and. op%const_w)
then
671 call accel_create_buffer(op%buff_weights, accel_mem_read_only, type_float, op%stencil%size)
672 call accel_create_buffer(op%buff_half_weights, accel_mem_read_only, type_float, op%stencil%size)
687 if (accel_is_enabled() .and. op%const_w)
then
688 call accel_write_buffer(op%buff_weights, op%stencil%size, op%w(:, 1))
689 call accel_write_buffer(op%buff_half_weights, op%stencil%size, -m_half*op%w(:, 1))
704 ii = op%rimap_inv(jj + 1) + maxval(op%ri(1:op%stencil%size, jj))
705 np_bc = max(np_bc, ii)
712#include "nl_operator_inc.F90"
715#include "complex.F90"
716#include "nl_operator_inc.F90"
subroutine, public accel_kernel_build(this, file_name, kernel_name, flags)
subroutine, public accel_release_buffer(this, async)
pure logical function, public accel_is_enabled()
integer, parameter, public accel_mem_read_only
integer pure function, public accel_max_workgroup_size()
This module implements batches of mesh functions.
Module implementing boundary conditions in Octopus.
real(real64), parameter, public m_zero
real(real64), parameter, public m_half
This module implements the index, used for the mesh points.
This module is intended to contain "only mathematical" functions and procedures.
This module defines the meshes, which are used in Octopus.
integer function, public mesh_local_index_from_coords(mesh, ix)
This function returns the local index of the point for a given vector of integer coordinates.
subroutine, public mesh_local_index_to_coords(mesh, ip, ix)
Given a local point index, this function returns the set of integer coordinates of the point.
subroutine, public messages_obsolete_variable(namespace, name, rep)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_experimental(name, namespace)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
This module handles the communicators for the various parallelization strategies.
This module defines non-local operators.
subroutine, public dnl_operator_operate_diag(op, fo)
subroutine, public nl_operator_init(op, label)
initialize an instance of a non-local operator by setting the label
integer, parameter op_map
integer, parameter op_max
integer, parameter op_vec
subroutine, public dnl_operator_operate_batch(op, fi, fo, ghost_update, profile, points, factor, async)
subroutine, public dnl_operator_operate(op, fi, fo, ghost_update, profile, points)
subroutine, public nl_operator_update_gpu_buffers(op)
subroutine, public nl_operator_global_init(namespace)
initialize global settings for non-local operators
subroutine, public nl_operator_copy(opo, opi)
subroutine, public nl_operator_output_weights(this)
subroutine, public nl_operator_end(op)
integer, parameter, public op_inner
subroutine, public znl_operator_operate(op, fi, fo, ghost_update, profile, points)
subroutine, public nl_operator_global_end()
integer, parameter, public op_outer
subroutine, public nl_operator_build(space, mesh, op, np, const_w, regenerate)
subroutine, public znl_operator_operate_batch(op, fi, fo, ghost_update, profile, points, factor, async)
integer, parameter op_nomap
integer pure function, public nl_operator_np_zero_bc(op)
subroutine, public znl_operator_operate_diag(op, fo)
integer pure function, public nl_operator_get_index(op, is, ip)
subroutine, public nl_operator_allocate_gpu_buffers(op)
integer, parameter op_min
This module contains interfaces for routines in operate.c.
Some general things and nomenclature:
This module defines stencils used in Octopus.
subroutine, public stencil_end(this)
subroutine, public stencil_copy(input, output)
type(type_t), public type_float
type(type_t), public type_integer
Describes mesh distribution to nodes.
index type for non-local operators
data type for non local operators