46 use json_module,
only : json_file, json_value, json_core
57 use,
intrinsic :: iso_c_binding
63 integer :: n_fields = 0
66 integer :: n_global_probes
68 integer :: n_probes_offset
70 real(kind=
rp),
allocatable :: xyz(:,:)
72 real(kind=
rp),
allocatable :: out_values(:,:)
73 type(c_ptr),
allocatable :: out_values_d(:)
74 real(kind=
rp),
allocatable :: out_vals_trsp(:,:)
76 integer :: n_local_probes
79 character(len=20),
allocatable :: which_fields(:)
81 integer,
allocatable :: n_local_probes_tot_offset(:)
82 integer,
allocatable :: n_local_probes_tot(:)
85 real(kind=
rp),
allocatable :: global_output_values(:,:)
93 procedure, pass(this) :: init_from_attributes => &
124 class(
probes_t),
intent(inout) :: this
125 type(json_file),
intent(inout) :: json
126 class(
case_t),
intent(inout),
target :: case
127 character(len=:),
allocatable :: output_file
128 character(len=:),
allocatable :: input_file
132 character(len=:),
allocatable :: point_type
133 type(json_value),
pointer :: json_point_list
134 type(json_file) :: json_point
135 type(json_core) :: core
136 integer :: idx, n_point_children
140 call this%init_base(json,
case)
143 call json%info(
'fields', n_children = this%n_fields)
144 call json_get(json,
'fields', this%which_fields)
145 call json_get(json,
'output_file', output_file)
147 call this%sampled_fields%init(this%n_fields)
148 do i = 1, this%n_fields
150 call this%sampled_fields%assign(i, &
155 this%n_local_probes = 0
156 this%n_global_probes = 0
159 if (json%valid_path(
'points_file'))
then
162 call json_get(json,
'points_file', input_file)
167 this%n_global_probes, input_file)
171 call json%get(
'points', json_point_list)
172 call json%info(
'points', n_children = n_point_children)
174 do idx = 1, n_point_children
178 select case (point_type)
181 call this%read_file(json_point)
183 call this%read_point(json_point)
185 call this%read_line(json_point)
187 call neko_error(
'Plane probes not implemented yet.')
189 call this%read_circle(json_point)
192 call this%read_point_zone(json_point,
case%fluid%dm_Xh)
195 call json_point%print()
198 call neko_error(
'Unknown region type ' // point_type)
202 call mpi_allreduce(this%n_local_probes, this%n_global_probes, 1, &
206 call this%init_from_attributes(
case%fluid%dm_Xh, output_file)
218 class(
probes_t),
intent(inout) :: this
219 type(json_file),
intent(inout) :: json
220 character(len=:),
allocatable :: input_file
221 real(kind=
rp),
dimension(:,:),
allocatable :: point_list
223 integer :: n_local, n_global
227 call json_get(json,
'file_name', input_file)
231 call this%add_points(point_list)
239 class(
probes_t),
intent(inout) :: this
240 type(json_file),
intent(inout) :: json
242 real(kind=
rp),
dimension(:,:),
allocatable :: point_list
243 real(kind=
rp),
dimension(:),
allocatable :: rp_list_reader
248 call json_get(json,
'coordinates', rp_list_reader)
251 if (.not. found)
then
255 if (mod(
size(rp_list_reader), 3) /= 0)
then
256 call neko_error(
'Invalid number of coordinates.')
260 allocate(point_list(3,
size(rp_list_reader)/3))
261 point_list = reshape(rp_list_reader, [3,
size(rp_list_reader)/3])
263 call this%add_points(point_list)
270 class(
probes_t),
intent(inout) :: this
271 type(json_file),
intent(inout) :: json
273 real(kind=
rp),
dimension(:,:),
allocatable :: point_list
274 real(kind=
rp),
dimension(:),
allocatable :: start,
end
275 real(kind=
rp),
dimension(3) :: direction
278 integer :: n_points, i
284 call json_get(json,
"amount", n_points)
287 if (
size(start) /= 3 .or.
size(
end) /= 3) then
288 call neko_error(
'Invalid start or end coordinates.')
292 allocate(point_list(3, n_points))
295 direction =
end - start
297 t =
real(i - 1, kind = rp) /
real(n_points - 1, kind = rp)
298 point_list(:, i) = start + direction * t
301 call this%add_points(point_list)
314 class(
probes_t),
intent(inout) :: this
315 type(json_file),
intent(inout) :: json
317 real(kind=rp),
dimension(:,:),
allocatable :: point_list
318 real(kind=rp),
dimension(:),
allocatable :: center, normal
319 real(kind=rp) :: radius
320 real(kind=rp) :: angle
321 integer :: n_points, i
322 character(len=:),
allocatable :: axis
324 real(kind=rp),
dimension(3) :: zero_line, cross_line, temp
328 if (pe_rank .ne. 0)
return
329 call json_get(json,
"center", center)
330 call json_get(json,
"normal", normal)
331 call json_get(json,
"radius", radius)
332 call json_get(json,
"amount", n_points)
333 call json_get(json,
"axis", axis)
336 if (
size(center) /= 3 .or.
size(normal) /= 3)
then
337 call neko_error(
'Invalid center or normal coordinates.')
339 if (axis /=
'x' .and. axis /=
'y' .and. axis /=
'z')
then
340 call neko_error(
'Invalid axis.')
342 if (radius <= 0)
then
343 call neko_error(
'Invalid radius.')
345 if (n_points <= 0)
then
346 call neko_error(
'Invalid number of points.')
350 normal = normal / norm2(normal)
353 if (axis .eq.
'x') zero_line = [1.0, 0.0, 0.0]
354 if (axis .eq.
'y') zero_line = [0.0, 1.0, 0.0]
355 if (axis .eq.
'z') zero_line = [0.0, 0.0, 1.0]
357 if (1.0_rp - dot_product(zero_line, normal) .le. 1e-6)
then
358 call neko_error(
'Invalid axis and normal.')
361 zero_line = zero_line - dot_product(zero_line, normal) * normal
362 zero_line = zero_line / norm2(zero_line)
364 cross_line(1) = normal(2) * zero_line(3) - normal(3) * zero_line(2)
365 cross_line(2) = normal(3) * zero_line(1) - normal(1) * zero_line(3)
366 cross_line(3) = normal(1) * zero_line(2) - normal(2) * zero_line(1)
369 allocate(point_list(3, n_points))
371 pi = 4.0_rp * atan(1.0_rp)
375 angle = 2.0_rp * pi *
real(i - 1, kind = rp) /
real(n_points, kind = rp)
376 temp = cos(angle) * zero_line + sin(angle) * cross_line
378 point_list(:, i) = center + radius * temp
381 call this%add_points(point_list)
390 class(
probes_t),
intent(inout) :: this
391 type(json_file),
intent(inout) :: json
392 type(dofmap_t),
intent(in) :: dof
394 real(kind=rp),
dimension(:,:),
allocatable :: point_list
395 character(len=:),
allocatable :: point_zone_name
396 class(point_zone_t),
pointer :: zone
397 integer :: i, idx, lx, nlindex(4)
398 real(kind=rp) :: x, y, z
401 if (pe_rank .ne. 0)
return
403 call json_get(json,
"name", point_zone_name)
404 zone => neko_point_zone_registry%get_point_zone(point_zone_name)
407 allocate(point_list(3, zone%size))
414 x = dof%x(nlindex(1), nlindex(2), nlindex(3), nlindex(4))
415 y = dof%y(nlindex(1), nlindex(2), nlindex(3), nlindex(4))
416 z = dof%z(nlindex(1), nlindex(2), nlindex(3), nlindex(4))
418 point_list(:, i) = [x, y, z]
421 call this%add_points(point_list)
431 class(
probes_t),
intent(inout) :: this
432 real(kind=rp),
dimension(:,:),
intent(in) :: new_points
434 real(kind=rp),
dimension(:,:),
allocatable :: temp
435 integer :: n_old, n_new
438 n_old = this%n_local_probes
439 n_new =
size(new_points, 2)
442 if (
allocated(this%xyz))
then
443 call move_alloc(this%xyz, temp)
447 allocate(this%xyz(3, n_old + n_new))
448 if (
allocated(temp))
then
449 this%xyz(:, 1:n_old) = temp
451 this%xyz(:, n_old+1:n_old+n_new) = new_points
453 this%n_local_probes = this%n_local_probes + n_new
463 class(
probes_t),
intent(inout) :: this
464 type(dofmap_t),
intent(in) :: dof
465 character(len=:),
allocatable,
intent(inout) :: output_file
466 character(len=1024) :: header_line
467 real(kind=rp),
allocatable :: global_output_coords(:,:)
469 type(matrix_t) :: mat_coords
472 call this%global_interp%init(dof)
475 call this%global_interp%find_points_and_redist(this%xyz, &
479 allocate(this%out_values(this%n_local_probes, this%n_fields))
480 allocate(this%out_values_d(this%n_fields))
481 allocate(this%out_vals_trsp(this%n_fields, this%n_local_probes))
483 if (neko_bcknd_device .eq. 1)
then
484 do i = 1, this%n_fields
485 this%out_values_d(i) = c_null_ptr
486 call device_map(this%out_values(:,i), this%out_values_d(i), &
492 this%fout = file_t(trim(output_file))
494 select type (ft => this%fout%file_type)
500 write(header_line,
'(I0,A,I0)') this%n_global_probes,
",", this%n_fields
501 do i = 1, this%n_fields
502 header_line = trim(header_line) //
"," // trim(this%which_fields(i))
504 call this%fout%set_header(header_line)
509 allocate(this%n_local_probes_tot(pe_size))
510 allocate(this%n_local_probes_tot_offset(pe_size))
511 call this%setup_offset()
512 if (pe_rank .eq. 0)
then
513 allocate(global_output_coords(3, this%n_global_probes))
514 call this%mat_out%init(this%n_global_probes, this%n_fields)
515 allocate(this%global_output_values(this%n_fields, &
516 this%n_global_probes))
517 call mat_coords%init(this%n_global_probes,3)
519 call mpi_gatherv(this%xyz, 3*this%n_local_probes, &
520 mpi_double_precision, global_output_coords, &
521 3*this%n_local_probes_tot, &
522 3*this%n_local_probes_tot_offset, &
523 mpi_double_precision, 0, neko_comm, ierr)
524 if (pe_rank .eq. 0)
then
525 call trsp(mat_coords%x, this%n_global_probes, &
526 global_output_coords, 3)
528 call this%fout%write(mat_coords)
531 call neko_error(
"Invalid data. Expected csv_file_t.")
538 class(
probes_t),
intent(inout) :: this
540 if (
allocated(this%xyz))
then
544 if (
allocated(this%out_values))
then
545 deallocate(this%out_values)
548 if (
allocated(this%out_vals_trsp))
then
549 deallocate(this%out_vals_trsp)
552 call this%sampled_fields%free()
554 if (
allocated(this%n_local_probes_tot))
then
555 deallocate(this%n_local_probes_tot)
558 if (
allocated(this%n_local_probes_tot_offset))
then
559 deallocate(this%n_local_probes_tot_offset)
562 if (
allocated(this%global_output_values))
then
563 deallocate(this%global_output_values)
566 call this%global_interp%free()
573 character(len=LOG_SIZE) :: log_buf
577 call neko_log%section(
'Probes')
578 write(log_buf,
'(A,I6)')
"Number of probes: ", this%n_global_probes
579 call neko_log%message(log_buf)
581 call neko_log%message(
"xyz-coordinates:", lvl = neko_log_debug)
582 do i = 1, this%n_local_probes
583 write(log_buf,
'("(",F10.6,",",F10.6,",",F10.6,")")') this%xyz(:,i)
584 call neko_log%message(log_buf, lvl = neko_log_debug)
588 write(log_buf,
'(A,I6)')
"Number of fields: ", this%n_fields
589 call neko_log%message(log_buf)
590 do i = 1, this%n_fields
591 write(log_buf,
'(A,I6,A,A)') &
592 "Field: ", i,
" ", trim(this%which_fields(i))
593 call neko_log%message(log_buf, lvl = neko_log_debug)
595 call neko_log%end_section()
596 call neko_log%newline()
604 character(len=LOG_SIZE) :: log_buf
607 do i = 1, this%n_local_probes
608 write (log_buf, *) pe_rank,
"/", this%global_interp%proc_owner(i), &
609 "/" , this%global_interp%el_owner(i), &
610 "/", this%global_interp%error_code(i)
611 call neko_log%message(log_buf)
612 write(log_buf,
'(A5,"(",F10.6,",",F10.6,",",F10.6,")")') &
613 "rst: ", this%global_interp%rst(:,i)
614 call neko_log%message(log_buf)
622 this%n_local_probes_tot = 0
623 this%n_local_probes_tot_offset = 0
624 this%n_probes_offset = 0
625 call mpi_gather(this%n_local_probes, 1, mpi_integer, &
626 this%n_local_probes_tot, 1, mpi_integer, &
629 call mpi_exscan(this%n_local_probes, this%n_probes_offset, 1, &
630 mpi_integer, mpi_sum, neko_comm, ierr)
631 call mpi_gather(this%n_probes_offset, 1, mpi_integer, &
632 this%n_local_probes_tot_offset, 1, mpi_integer, &
644 class(
probes_t),
intent(inout) :: this
645 real(kind=rp),
intent(in) :: t
646 integer,
intent(in) :: tstep
650 do i = 1, this%n_fields
651 call this%global_interp%evaluate(this%out_values(:,i), &
652 this%sampled_fields%items(i)%ptr%x)
655 if (neko_bcknd_device .eq. 1)
then
656 do i = 1, this%n_fields
657 call device_memcpy(this%out_values(:,i), this%out_values_d(i), &
658 this%n_local_probes, device_to_host, sync = .true.)
662 if (this%output_controller%check(t, tstep))
then
665 if (this%seq_io)
then
666 call trsp(this%out_vals_trsp, this%n_fields, &
667 this%out_values, this%n_local_probes)
668 call mpi_gatherv(this%out_vals_trsp, &
669 this%n_fields*this%n_local_probes, &
670 mpi_double_precision, this%global_output_values, &
671 this%n_fields*this%n_local_probes_tot, &
672 this%n_fields*this%n_local_probes_tot_offset, &
673 mpi_double_precision, 0, neko_comm, ierr)
674 if (pe_rank .eq. 0)
then
675 call trsp(this%mat_out%x, this%n_global_probes, &
676 this%global_output_values, this%n_fields)
677 call this%fout%write(this%mat_out, t)
680 call neko_error(
'probes sim comp, parallel io need implementation')
684 call this%output_controller%register_execution()
693 class(
probes_t),
intent(inout) :: this
694 character(len=:),
allocatable :: points_file
695 real(kind=rp),
allocatable :: xyz(:,:)
696 integer,
intent(inout) :: n_local_probes, n_global_probes
699 type(file_t) :: file_in
701 file_in = file_t(trim(points_file))
703 select type (ft => file_in%file_type)
708 call neko_error(
"Invalid data. Expected csv_file_t.")
712 call file_free(file_in)
722 type(csv_file_t),
intent(inout) :: f
723 real(kind=rp),
allocatable :: xyz(:,:)
724 integer,
intent(inout) :: n_local_probes, n_global_probes
725 type(matrix_t) :: mat_in, mat_in2
728 n_lines = f%count_lines()
731 n_global_probes = n_lines
734 if (pe_rank .eq. 0)
then
735 n_local_probes = n_global_probes
736 allocate(xyz(3, n_local_probes))
737 call mat_in%init(n_global_probes,3)
738 call mat_in2%init(3, n_global_probes)
740 call trsp(xyz, 3, mat_in%x, n_global_probes)
743 allocate(xyz(3, n_local_probes))
__device__ void nonlinear_index(const int idx, const int lx, int *index)
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Retrieves a parameter by name or throws an error.
Defines a simulation case.
type(mpi_comm) neko_comm
MPI communicator.
File format for .csv files, used for any read/write operations involving floating point data.
Device abstraction, common interface for various accelerators.
Defines a mapping of the degrees of freedom.
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Module for file I/O operations.
subroutine file_free(this)
File operation destructor.
Implements global_interpolation given a dofmap.
Utilities for retrieving parameters from the case files.
integer, parameter, public neko_log_debug
Debug log level.
type(log_t), public neko_log
Global log stream.
integer, parameter, public log_size
integer, parameter, public rp
Global precision used in computations.
type(point_zone_registry_t), target, public neko_point_zone_registry
Global point_zone registry.
subroutine read_point_zone(this, json, dof)
Construct a list of points from a point zone.
subroutine read_circle(this, json)
Construct a list of points from a circle.
subroutine add_points(this, new_points)
Append a new list of points to the exsiting list.
subroutine probes_init_from_attributes(this, dof, output_file)
Initialize without json things.
subroutine read_line(this, json)
Construct a list of points from a line.
subroutine read_probe_locations(this, xyz, n_local_probes, n_global_probes, points_file)
Initialize the physical coordinates from a csv input file.
subroutine read_point(this, json)
Read a list of points from the json file.
subroutine probes_show(this)
Print current probe status, with number of probes and coordinates.
subroutine probes_free(this)
Destructor.
subroutine probes_debug(this)
Show the status of processor/element owner and error code for each point.
subroutine probes_setup_offset(this)
Setup offset for rank 0.
subroutine probes_evaluate_and_write(this, t, tstep)
Interpolate each probe from its r,s,t coordinates.
subroutine read_file(this, json)
Read a list of points from a csv file.
subroutine read_xyz_from_csv(xyz, n_local_probes, n_global_probes, f)
Read and initialize the number of probes from a csv input file.
subroutine probes_init_from_json(this, json, case)
Constructor from json.
Simulation components are objects that encapsulate functionality that can be fit to a particular comp...
subroutine, public trsp(a, lda, b, ldb)
Transpose of a rectangular tensor .
field_list_t, To be able to group fields together
A wrapper around a polymorphic generic_file_t that handles its init. This is essentially a factory fo...
Implements global interpolation for arbitrary points in the domain.
Base abstract type for point zones.
Base abstract class for simulation components.