Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
global_interpolation.f90
Go to the documentation of this file.
1! Copyright (c) 2020-2025, 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!
36 use num_types, only: rp, dp, xp
38 use space, only: space_t
39 use stack, only: stack_i4_t
40 use dofmap, only: dofmap_t
41 use logger, only: neko_log, log_size
42 use utils, only: neko_error
52 use el_finder, only: el_finder_t
53 use pe_finder, only: pe_finder_t
54 use comm, only: neko_comm
55 use mpi_f08, only: mpi_sum, mpi_comm, mpi_comm_rank, &
56 mpi_comm_size, mpi_wtime, mpi_allreduce, mpi_in_place, mpi_integer, &
57 mpi_min, mpi_barrier, mpi_reduce_scatter_block, mpi_alltoall, &
58 mpi_isend, mpi_irecv
60 use vector, only: vector_t
61 use matrix, only: matrix_t
62 use math, only: copy, neko_eps
63 use structs, only : array_ptr_t
64 use, intrinsic :: iso_c_binding, only: c_ptr, c_null_ptr, c_associated
65 implicit none
66 private
67
68 integer, public, parameter :: glob_map_size = 4096
69
71 type, public :: global_interpolation_t
73 type(vector_t) :: x
75 type(vector_t) :: y
77 type(vector_t) :: z
79 integer :: gdim
81 integer :: nelv
83 integer :: glb_nelv
85 type(mpi_comm) :: comm
87 integer :: pe_rank
89 integer :: pe_size
91 type(space_t) :: xh
94 integer :: n_points
96 real(kind=rp), allocatable :: xyz(:,:)
98 real(kind=rp), allocatable :: rst(:,:)
100 integer, allocatable :: pe_owner(:)
102 type(stack_i4_t), allocatable :: points_at_pe(:)
105 integer, allocatable :: el_owner0(:)
106 type(c_ptr) :: el_owner0_d = c_null_ptr
107
109 integer :: n_points_local
110 integer, allocatable :: el_owner0_local(:)
111 type(c_ptr) :: el_owner0_local_d = c_null_ptr
112 real(kind=rp), allocatable :: rst_local(:,:)
113 real(kind=rp), allocatable :: xyz_local(:,:)
115 type(local_interpolator_t) :: local_interp
118 logical :: all_points_local = .false.
120 real(kind=rp) :: tol = neko_eps*1e3_rp
122 real(kind=rp) :: padding = 1e-2_rp
123
127 integer, allocatable :: n_points_pe(:)
128 integer, allocatable :: n_points_offset_pe(:)
131 integer, allocatable :: n_points_pe_local(:)
132 integer, allocatable :: n_points_offset_pe_local(:)
135 class(pe_finder_t), allocatable :: pe_finder
137 class(el_finder_t), allocatable :: el_finder
139 type(legendre_rst_finder_t) :: rst_finder
143 type(vector_t) :: temp_local, temp
144
145 contains
147 procedure, pass(this) :: init_xyz => global_interpolation_init_xyz
149 procedure, pass(this) :: init_dof => global_interpolation_init_dof
151 procedure, pass(this) :: free => global_interpolation_free
153 procedure, pass(this) :: free_points => global_interpolation_free_points
154 procedure, pass(this) :: free_points_local => global_interpolation_free_points_local
155 procedure, pass(this) :: find_points_and_redist => &
160 procedure, pass(this) :: find_points_coords => &
162 procedure, pass(this) :: find_points_coords1d => &
165 procedure, pass(this) :: check_points => &
167 procedure, pass(this) :: find_points_xyz => global_interpolation_find_xyz
168 generic :: find_points => find_points_xyz, find_points_coords, find_points_coords1d
170 procedure, pass(this) :: evaluate => global_interpolation_evaluate
171
173 generic :: init => init_dof, init_xyz
174
176
177contains
178
182 subroutine global_interpolation_init_dof(this, dof, comm, tol, pad)
183 class(global_interpolation_t), target, intent(inout) :: this
184 type(dofmap_t) :: dof
185 type(mpi_comm), optional, intent(in) :: comm
186 real(kind=rp), optional :: tol
187 real(kind=rp), optional :: pad
188 real(kind=rp) :: padding, tolerance
189
190 if (present(pad)) then
191 padding = pad
192 else
193 padding = 1e-2 ! 1% padding of the bounding boxes
194 end if
195
196 if (present(tol)) then
197 tolerance = tol
198 else
199 tolerance = neko_eps*1e3_xp ! 1% padding of the bounding boxes
200 end if
201 ! NOTE: Passing dof%x(:,1,1,1), etc in init_xyz passes down the entire
202 ! dof%x array and not a slice. It is done this way for
203 ! to get the right dimension (see global_interpolation_init_xyz).
204 call this%init_xyz(dof%x(:,1,1,1), dof%y(:,1,1,1), dof%z(:,1,1,1), &
205 dof%msh%gdim, dof%msh%nelv, dof%Xh, comm,tol = tolerance, pad=padding)
206
207 end subroutine global_interpolation_init_dof
208
219 subroutine global_interpolation_init_xyz(this, x, y, z, gdim, nelv, Xh, &
220 comm, tol, pad)
221 class(global_interpolation_t), target, intent(inout) :: this
222 real(kind=rp), intent(in) :: x(:)
223 real(kind=rp), intent(in) :: y(:)
224 real(kind=rp), intent(in) :: z(:)
225 integer, intent(in) :: gdim
226 integer, intent(in) :: nelv
227 type(mpi_comm), intent(in), optional :: comm
228 type(space_t), intent(in) :: Xh
229 real(kind=rp), intent(in), optional :: tol
230 real(kind=rp), intent(in), optional :: pad
231 integer :: lx, ly, lz, ierr, i, n
232 character(len=8000) :: log_buf
233 real(kind=dp) :: padding
234 real(kind=rp) :: time1, time_start
235 character(len=255) :: mode_str
236 integer :: boxdim, envvar_len
237
238 call this%free()
239
240 if (present(comm)) then
241 this%comm = comm
242 else
243 this%comm = neko_comm
244 end if
245
246 if (present(pad)) then
247 padding = pad
248 this%padding = pad
249 else
250 padding = 1e-2 ! 1% padding of the bounding boxes
251 this%padding = 1e-2
252 end if
253
254 time_start = mpi_wtime()
255 call mpi_barrier(this%comm)
256
257 call mpi_comm_rank(this%comm, this%pe_rank, ierr)
258 call mpi_comm_size(this%comm, this%pe_size, ierr)
259
260 this%gdim = gdim
261 this%nelv = nelv
262 if (present(tol)) then
263 this%tol = tol
264 else
265 this%tol = neko_eps*1e3_xp
266 end if
267
268 call mpi_allreduce(nelv, this%glb_nelv, 1, mpi_integer, &
269 mpi_sum, this%comm, ierr)
270 lx = xh%lx
271 ly = xh%ly
272 lz = xh%lz
273 n = nelv * lx*ly*lz
274 call this%x%init(n)
275 call this%y%init(n)
276 call this%z%init(n)
277 call copy(this%x%x, x, n)
278 call this%x%copy_from(host_to_device,.false.)
279 call copy(this%y%x, y, n)
280 call this%y%copy_from(host_to_device,.false.)
281 call copy(this%z%x, z, n)
282 call this%z%copy_from(host_to_device,.false.)
283 call this%Xh%init(xh%t, lx, ly, lz)
284
285
286 call neko_log%section('Global Interpolation')
287
288 call neko_log%message('Initializing global interpolation')
289 call get_environment_variable("NEKO_GLOBAL_INTERP_EL_FINDER", &
290 mode_str, envvar_len)
291
292 if (envvar_len .gt. 0) then
293 if (mode_str(1:envvar_len) == 'AABB') then
294 allocate(aabb_el_finder_t :: this%el_finder)
295 end if
296 end if
297 call get_environment_variable("NEKO_GLOBAL_INTERP_PE_FINDER", &
298 mode_str, envvar_len)
299
300 if (envvar_len .gt. 0) then
301 if (mode_str(1:envvar_len) == 'AABB') then
302 allocate(aabb_pe_finder_t :: this%pe_finder)
303 end if
304 end if
305
306 if (.not. allocated(this%el_finder)) then
307 allocate(cartesian_el_finder_t :: this%el_finder)
308 end if
309 if (.not. allocated(this%pe_finder)) then
310 allocate(cartesian_pe_finder_t :: this%pe_finder)
311 end if
312 select type(el_find => this%el_finder)
313 type is (aabb_el_finder_t)
314 call neko_log%message('Using AABB element finder')
315 call el_find%init(x, y, z, nelv, xh, padding)
316 type is (cartesian_el_finder_t)
317 call neko_log%message('Using Cartesian element finder')
318 boxdim = max(lx*int(real(nelv,xp)**(1.0_xp/3.0_xp)),2)
319 boxdim = min(boxdim, 300)
320 call el_find%init(x, y, z, nelv, xh, boxdim, padding)
321 class default
322 call neko_error('Unknown element finder type')
323 end select
324
325 select type(pe_find => this%pe_finder)
326 type is (aabb_pe_finder_t)
327 call neko_log%message('Using AABB PE finder')
328 call pe_find%init(this%x%x, this%y%x, this%z%x, &
329 nelv, xh, this%comm, padding)
330 type is (cartesian_pe_finder_t)
331 call neko_log%message('Using Cartesian PE finder')
332 boxdim = lx*int(real(this%glb_nelv,xp)**(1.0_xp/3.0_xp))
333 boxdim = max(boxdim,32)
334 boxdim = min(boxdim, &
335 int(8.0_xp*(30000.0_xp*this%pe_size)**(1.0_xp/3.0_xp)))
336 call pe_find%init(this%x%x, this%y%x, this%z%x, &
337 nelv, xh, this%comm, boxdim, padding)
338 class default
339 call neko_error('Unknown PE finder type')
340 end select
341
342 call this%rst_finder%init(this%x%x, this%y%x, this%z%x, nelv, xh, this%tol)
343 if (allocated(this%n_points_pe)) deallocate(this%n_points_pe)
344 if (allocated(this%n_points_pe_local)) deallocate(this%n_points_pe_local)
345 if (allocated(this%n_points_offset_pe_local)) &
346 deallocate(this%n_points_offset_pe_local)
347 if (allocated(this%n_points_offset_pe)) deallocate(this%n_points_offset_pe)
348 allocate(this%n_points_pe(0:(this%pe_size-1)))
349 allocate(this%n_points_offset_pe(0:(this%pe_size-1)))
350 allocate(this%n_points_pe_local(0:(this%pe_size-1)))
351 allocate(this%n_points_offset_pe_local(0:(this%pe_size-1)))
352 allocate(this%points_at_pe(0:(this%pe_size-1)))
353 do i = 0, this%pe_size-1
354 call this%points_at_pe(i)%init()
355 end do
356 call mpi_barrier(this%comm)
357 time1 = mpi_wtime()
358 write(log_buf, '(A,E15.7)') &
359 'Global interpolation initialized (s):', time1-time_start
360 call neko_log%message(log_buf)
361 call neko_log%end_section()
362 end subroutine global_interpolation_init_xyz
363
364
367 class(global_interpolation_t), target, intent(inout) :: this
368 integer :: i
369
370 call this%x%free()
371 call this%y%free()
372 call this%z%free()
373 call this%Xh%free()
374
375 this%nelv = 0
376 this%gdim = 0
377
378 call this%free_points()
379 call this%free_points_local()
380 call this%local_interp%free()
381 if (allocated(this%el_finder)) then
382 call this%el_finder%free()
383 deallocate(this%el_finder)
384 end if
385 if (allocated(this%pe_finder)) then
386 call this%pe_finder%free()
387 deallocate(this%pe_finder)
388 end if
389 call this%rst_finder%free()
390
391 call this%temp_local%free()
392 call this%temp%free()
393 if (allocated(this%points_at_pe)) then
394 do i = 0, this%pe_size-1
395 call this%points_at_pe(i)%free()
396 end do
397 deallocate(this%points_at_pe)
398 end if
399 if (allocated(this%n_points_pe)) deallocate(this%n_points_pe)
400 if (allocated(this%n_points_pe_local)) deallocate(this%n_points_pe_local)
401 if (allocated(this%n_points_offset_pe_local)) &
402 deallocate(this%n_points_offset_pe_local)
403 if (allocated(this%n_points_offset_pe)) deallocate(this%n_points_offset_pe)
404
405
406
407 end subroutine global_interpolation_free
408
411 class(global_interpolation_t), target, intent(inout) :: this
412
413 this%n_points = 0
414 this%all_points_local = .false.
415
416 if (allocated(this%xyz)) deallocate(this%xyz)
417 if (allocated(this%rst)) deallocate(this%rst)
418 if (allocated(this%pe_owner)) deallocate(this%pe_owner)
419 if (allocated(this%el_owner0)) deallocate(this%el_owner0)
420
421 if (c_associated(this%el_owner0_d)) then
422 call device_free(this%el_owner0_d)
423 end if
424
425 call this%glb_intrp_comm%free()
426
427
429
430
432 class(global_interpolation_t), target, intent(inout) :: this
433
434 this%n_points_local = 0
435 this%all_points_local = .false.
436
437 if (allocated(this%xyz_local)) deallocate(this%xyz_local)
438 if (allocated(this%rst_local)) deallocate(this%rst_local)
439 if (allocated(this%el_owner0_local)) deallocate(this%el_owner0_local)
440
441 if (c_associated(this%el_owner0_local_d)) then
442 call device_free(this%el_owner0_local_d)
443 end if
444
446
447
450 class(global_interpolation_t), target, intent(inout) :: this
451 character(len=8000) :: log_buf
452 type(vector_t) :: x_t
453 type(vector_t) :: y_t
454 type(vector_t) :: z_t
455 type(matrix_t) :: rst_local_cand
456 type(vector_t) :: resx
457 type(vector_t) :: resy
458 type(vector_t) :: resz
459 type(c_ptr) :: el_cands_d
460 type(matrix_t) :: res
461 integer :: i, j, stupid_intent
462 type(stack_i4_t), target :: all_el_candidates
463 integer, allocatable :: n_el_cands(:)
464 integer, contiguous, pointer :: el_cands(:), point_ids(:), send_recv(:)
465 real(kind=rp), allocatable :: res_results(:,:)
466 real(kind=rp), allocatable :: rst_results(:,:)
467 integer, allocatable :: el_owner_results(:)
468 integer :: ierr, ii, n_point_cand, n_glb_point_cand, point_id, rank
469 real(kind=rp) :: time1, time2, time_start
470 !Temp stuff for glb_intrp_comm
471 type(stack_i4_t) :: send_pe, recv_pe
472 type(glb_intrp_comm_t) :: glb_intrp_find, glb_intrp_find_back
473 type(stack_i4_t) :: send_pe_find, recv_pe_find
474
475 el_cands_d = c_null_ptr
476 call neko_log%section('Global Interpolation')
477 call glb_intrp_find%init_dofs(this%pe_size)
478 call send_pe_find%init()
479 call recv_pe_find%init()
480 call mpi_barrier(this%comm)
481 time_start = mpi_wtime()
482 write(log_buf,'(A)') 'Global interpolation, finding points'
483 call neko_log%message(log_buf)
484 ! Find pe candidates that the points i want may be at
485 ! Add number to n_points_pe_local
486 !Working arrays
487 this%n_points_pe = 0
488 call this%pe_finder%find_batch(this%xyz, this%n_points, &
489 this%points_at_pe, this%n_points_pe)
490 call mpi_barrier(this%comm)
491 time1 = mpi_wtime()
492 write(log_buf, '(A,E15.7)') &
493 'Found PE candidates time since start of findpts (s):', time1-time_start
494 call neko_log%message(log_buf)
495
496 !Send number of points I want to candidates
497 ! n_points_local -> how many points might be at this rank
498 ! n_points_pe_local -> how many points local on this rank that other pes might want
499 this%n_points_pe_local = 0
500 this%n_points_local = 0
501 call mpi_reduce_scatter_block(this%n_points_pe, this%n_points_local, &
502 1, mpi_integer, mpi_sum, this%comm, ierr)
503 call mpi_alltoall(this%n_points_pe, 1, mpi_integer,&
504 this%n_points_pe_local, 1, mpi_integer, this%comm, ierr)
505
506 !Set up offset arrays
507 this%n_points_offset_pe_local(0) = 0
508 this%n_points_offset_pe(0) = 0
509 do i = 1, (this%pe_size - 1)
510 this%n_points_offset_pe_local(i) = this%n_points_pe_local(i-1)&
511 + this%n_points_offset_pe_local(i-1)
512 this%n_points_offset_pe(i) = this%n_points_pe(i-1)&
513 + this%n_points_offset_pe(i-1)
514 end do
515 do i = 0, (this%pe_size-1)
516 if (this%n_points_pe(i) .gt. 0) then
517 call send_pe_find%push(i)
518 point_ids => this%points_at_pe(i)%array()
519 do j = 1, this%n_points_pe(i)
520 call glb_intrp_find%send_dof(i)%push(3*(point_ids(j)-1)+1)
521 call glb_intrp_find%send_dof(i)%push(3*(point_ids(j)-1)+2)
522 call glb_intrp_find%send_dof(i)%push(3*(point_ids(j)-1)+3)
523 end do
524 end if
525 if (this%n_points_pe_local(i) .gt. 0) then
526 call recv_pe_find%push(i)
527 do j = 1, this%n_points_pe_local(i)
528 call glb_intrp_find%recv_dof(i)%push(3*(j+this%n_points_offset_pe_local(i)-1)+1)
529 call glb_intrp_find%recv_dof(i)%push(3*(j+this%n_points_offset_pe_local(i)-1)+2)
530 call glb_intrp_find%recv_dof(i)%push(3*(j+this%n_points_offset_pe_local(i)-1)+3)
531 end do
532 end if
533 end do
534
535
536
537 call glb_intrp_find%init(send_pe_find, recv_pe_find, this%comm)
538
539 call glb_intrp_find_back%init_dofs(this%pe_size)
540 ii = 0
541 do i = 0, (this%pe_size-1)
542 send_recv => glb_intrp_find%recv_dof(i)%array()
543 do j = 1, glb_intrp_find%recv_dof(i)%size()
544 call glb_intrp_find_back%send_dof(i)%push(send_recv(j))
545 end do
546 send_recv => glb_intrp_find%send_dof(i)%array()
547 do j = 1, glb_intrp_find%send_dof(i)%size()
548 ii = ii + 1
549 call glb_intrp_find_back%recv_dof(i)%push(ii)
550 end do
551 end do
552
553 call glb_intrp_find_back%init(recv_pe_find, send_pe_find, this%comm)
554
555
556 if (allocated(this%xyz_local)) then
557 deallocate(this%xyz_local)
558 end if
559 allocate(this%xyz_local(3, this%n_points_local))
560 call glb_intrp_find%sendrecv(this%xyz, this%xyz_local, this%n_points*3, &
561 this%n_points_local*3)
562
563 call mpi_barrier(this%comm)
564 time1 = mpi_wtime()
565 write(log_buf, '(A,E15.7)') &
566 'Sent to points to PE candidates, time since start of find_points (s):', time1-time_start
567 call neko_log%message(log_buf)
568 !Okay, now we need to find the rst...
569 call all_el_candidates%init()
570
571 if (allocated(n_el_cands)) then
572 deallocate(n_el_cands)
573 end if
574
575 allocate(n_el_cands(this%n_points_local))
577 call this%el_finder%find_batch(this%xyz_local, this%n_points_local, &
578 all_el_candidates, n_el_cands)
579
580 n_point_cand = all_el_candidates%size()
581 if (n_point_cand .gt. 1e8) then
582 print *,'Warning, many point candidates on rank', this%pe_rank,'cands:', n_point_cand, &
583 'Consider increasing number of ranks'
584 end if
585 call x_t%init(n_point_cand)
586 call y_t%init(n_point_cand)
587 call z_t%init(n_point_cand)
588 ii = 0
590 do i = 1 , this%n_points_local
591 do j = 1, n_el_cands(i)
592 ii = ii + 1
593 x_t%x(ii) = this%xyz_local(1,i)
594 y_t%x(ii) = this%xyz_local(2,i)
595 z_t%x(ii) = this%xyz_local(3,i)
596 end do
597 end do
598
599 call mpi_barrier(this%comm)
600 time1 = mpi_wtime()
601 write(log_buf, '(A,E15.7)') &
602 'Element candidates found, now time for finding rst,time since start of find_points (s):', time1-time_start
603 call neko_log%message(log_buf)
604 call rst_local_cand%init(3,n_point_cand)
605 call resx%init(n_point_cand)
606 call resy%init(n_point_cand)
607 call resz%init(n_point_cand)
608
609 ! Find rst within all element candidates for target xyz (x_t, y_t, z_t)
610 call mpi_barrier(this%comm)
611 time1 = mpi_wtime()
612 el_cands => all_el_candidates%array()
613 if (neko_bcknd_device .eq. 1) then
614 call x_t%copy_from(host_to_device,.false.)
615
616 call y_t%copy_from(host_to_device,.false.)
617
618 call z_t%copy_from(host_to_device,.false.)
619 call device_map(el_cands, el_cands_d,n_point_cand)
620 call device_memcpy(el_cands, el_cands_d,n_point_cand, &
621 host_to_device, .true.)
622 end if
623
624 call this%rst_finder%find(rst_local_cand, &
625 x_t, y_t, z_t, &
626 el_cands, n_point_cand, &
627 resx, resy, resz)
628 if (neko_bcknd_device .eq. 1) then
629 call rst_local_cand%copy_from(device_to_host,.false.)
630 call resx%copy_from(device_to_host,.false.)
631 call resy%copy_from(device_to_host,.false.)
632 call resz%copy_from(device_to_host,.true.)
633 call device_deassociate(el_cands)
634 call device_free(el_cands_d)
635 end if
636 call mpi_barrier(this%comm)
637
638 time2 = mpi_wtime()
639 write(log_buf, '(A,E15.7)') &
640 'Found rst with Newton iteration, time (s):', time2-time1
641 call neko_log%message(log_buf)
642
643 write(log_buf,'(A,E15.7)') &
644 'Tolerance: ', this%tol
645 call neko_log%message(log_buf)
646 write(log_buf,'(A)') &
647 'Checking validity of points and choosing best candidates.'
648 call neko_log%message(log_buf)
649 call mpi_barrier(this%comm,ierr)
650
651 if (allocated(this%rst_local)) deallocate(this%rst_local)
652 if (allocated(this%el_owner0_local)) deallocate(this%el_owner0_local)
653 allocate(this%rst_local(3,this%n_points_local))
654 allocate(this%el_owner0_local(this%n_points_local))
655 ! Choose the best candidate at this rank
656 ii = 0
657 do i = 1 , this%n_points_local
658 this%xyz_local(1,i) = 10.0
659 this%xyz_local(2,i) = 10.0
660 this%xyz_local(3,i) = 10.0
661 this%rst_local(1,i) = 10.0
662 this%rst_local(2,i) = 10.0
663 this%rst_local(3,i) = 10.0
664 this%el_owner0_local(i) = -1
665 do j = 1, n_el_cands(i)
666 ii = ii + 1
667 if (rst_cmp(this%rst_local(:,i), rst_local_cand%x(:,ii),&
668 this%xyz_local(:,i), (/resx%x(ii),resy%x(ii),resz%x(ii)/), this%padding)) then
669 this%rst_local(1,i) = rst_local_cand%x(1,ii)
670 this%rst_local(2,i) = rst_local_cand%x(2,ii)
671
672 this%rst_local(3,i) = rst_local_cand%x(3,ii)
673 this%xyz_local(1,i) = resx%x(ii)
674 this%xyz_local(2,i) = resy%x(ii)
675 this%xyz_local(3,i) = resz%x(ii)
676 this%el_owner0_local(i) = el_cands(ii)
677 end if
678 ! if (this%pe_rank .eq. 0) print *,i, this%rst_local(:,i), &
679 ! this%xyz_local(:,i), this%el_owner0_local(i)
680 end do
681 end do
682 call res%init(3,this%n_points)
683 n_glb_point_cand = sum(this%n_points_pe)
684 if (allocated(rst_results)) deallocate(rst_results)
685 if (allocated(res_results)) deallocate(res_results)
686 if (allocated(el_owner_results)) deallocate(el_owner_results)
687 allocate(rst_results(3,n_glb_point_cand))
688 allocate(res_results(3,n_glb_point_cand))
689 allocate(el_owner_results(n_glb_point_cand))
690 res = 1e2_rp
691 this%rst = 1e2
692 this%pe_owner = -1
693 this%el_owner0 = -1
694 call glb_intrp_find_back%sendrecv(this%xyz_local, res_results, &
695 this%n_points_local*3, n_glb_point_cand*3)
696 call glb_intrp_find_back%sendrecv(this%rst_local, rst_results, &
697 this%n_points_local*3, n_glb_point_cand*3)
698 do i = 1, size(glb_intrp_find_back%send_pe)
699 rank = glb_intrp_find_back%send_pe(i)
700 call mpi_isend(this%el_owner0_local(this%n_points_offset_pe_local(rank)+1),&
701 this%n_points_pe_local(rank), &
702 mpi_integer, rank, 0, &
703 this%comm, glb_intrp_find_back%send_buf(i)%request, ierr)
704 glb_intrp_find_back%send_buf(i)%flag = .false.
705 end do
706 do i = 1, size(glb_intrp_find_back%recv_pe)
707 rank = glb_intrp_find_back%recv_pe(i)
708 call mpi_irecv(el_owner_results(this%n_points_offset_pe(rank)+1),&
709 this%n_points_pe(rank), &
710 mpi_integer, rank, 0, &
711 this%comm, glb_intrp_find_back%recv_buf(i)%request, ierr)
712 glb_intrp_find_back%recv_buf(i)%flag = .false.
713 end do
714 call glb_intrp_find_back%nbwait_no_op()
715 ii = 0
716 do i = 1, size(glb_intrp_find_back%recv_pe)
717 point_ids => this%points_at_pe(glb_intrp_find_back%recv_pe(i))%array()
718 do j = 1, this%n_points_pe(glb_intrp_find_back%recv_pe(i))
719 point_id = point_ids(j)
720 ii = ii + 1
721 if (rst_cmp(this%rst(:,point_id), rst_results(:,ii), &
722 res%x(:,point_id), res_results(:,ii), this%padding) .or. &
723 this%pe_owner(point_ids(j)) .eq. -1 ) then
724 this%rst(:,point_ids(j)) = rst_results(:,ii)
725 res%x(:,point_ids(j)) = res_results(:,ii)
726 this%pe_owner(point_ids(j)) = glb_intrp_find_back%recv_pe(i)
727 this%el_owner0(point_ids(j)) = el_owner_results(ii)
728 end if
729 ! if (this%pe_rank .eq. 0) print *,point_id, &
730 !this%rst(:,point_ids(j)),res%x(:,point_ids(j)), this%el_owner0(point_ids(j))
731 end do
732 end do
733
734 !OK, now I know the correct rst values of the points I want We now
735 !send the correct rsts to the correct rank (so a point only
736 !belongs to one rank)
737 do i = 0, this%pe_size-1
738 call this%points_at_pe(i)%clear()
739 this%n_points_pe(i) = 0
740 end do
741
742 do i = 1, this%n_points
743 stupid_intent = i
744 if (this%pe_owner(i) .eq. -1 .or. this%el_owner0(i) .eq. -1) then
745 print *, 'No owning rank found for',&
746 ' point ', stupid_intent, ' with coords', this%xyz(:,i), &
747 ' Interpolation will always yield 0.0. Try increase padding.'
748 else
749 call this%points_at_pe(this%pe_owner(i))%push(stupid_intent)
750 this%n_points_pe(this%pe_owner(i)) = this%n_points_pe(this%pe_owner(i)) + 1
751 end if
752 end do
753 call mpi_reduce_scatter_block(this%n_points_pe, this%n_points_local, 1, mpi_integer, &
754 mpi_sum, this%comm, ierr)
755 call mpi_alltoall(this%n_points_pe, 1, mpi_integer,&
756 this%n_points_pe_local, 1, mpi_integer, this%comm, ierr)
757 this%n_points_offset_pe_local(0) = 0
758 this%n_points_offset_pe(0) = 0
759 do i = 1, (this%pe_size - 1)
760 this%n_points_offset_pe_local(i) = this%n_points_pe_local(i-1)&
761 + this%n_points_offset_pe_local(i-1)
762 this%n_points_offset_pe(i) = this%n_points_pe(i-1)&
763 + this%n_points_offset_pe(i-1)
764 end do
765 call send_pe_find%free()
766 call recv_pe_find%free()
767 call glb_intrp_find%free()
768 call send_pe_find%init()
769 call recv_pe_find%init()
770 call glb_intrp_find%init_dofs(this%pe_size)
771 !setup comm to send xyz and rst to chosen ranks
772 do i = 0, (this%pe_size-1)
773 if (this%n_points_pe(i) .gt. 0) then
774 call send_pe_find%push(i)
775 point_ids => this%points_at_pe(i)%array()
776 do j = 1, this%n_points_pe(i)
777 call glb_intrp_find%send_dof(i)%push(3*(point_ids(j)-1)+1)
778 call glb_intrp_find%send_dof(i)%push(3*(point_ids(j)-1)+2)
779 call glb_intrp_find%send_dof(i)%push(3*(point_ids(j)-1)+3)
780 end do
781 end if
782 if (this%n_points_pe_local(i) .gt. 0) then
783 call recv_pe_find%push(i)
784 do j = 1, this%n_points_pe_local(i)
785 call glb_intrp_find%recv_dof(i)%push(3*(j+this%n_points_offset_pe_local(i)-1)+1)
786 call glb_intrp_find%recv_dof(i)%push(3*(j+this%n_points_offset_pe_local(i)-1)+2)
787 call glb_intrp_find%recv_dof(i)%push(3*(j+this%n_points_offset_pe_local(i)-1)+3)
788 end do
789 end if
790 end do
791
792
793 call glb_intrp_find%init(send_pe_find, recv_pe_find, this%comm)
794 call glb_intrp_find%sendrecv(this%xyz, this%xyz_local, this%n_points*3, &
795 this%n_points_local*3)
796 call glb_intrp_find%sendrecv(this%rst, this%rst_local, this%n_points*3, &
797 this%n_points_local*3)
798 ii = 0
799 do i = 1, size(glb_intrp_find%send_pe)
800 rank = glb_intrp_find%send_pe(i)
801 point_ids => this%points_at_pe(rank)%array()
802 do j = 1, this%n_points_pe(rank)
803 ii = ii + 1
804 el_owner_results(ii) = this%el_owner0(point_ids(j))
805 end do
806 call mpi_isend(el_owner_results(this%n_points_offset_pe(rank)+1),&
807 this%n_points_pe(rank), &
808 mpi_integer, rank, 0, &
809 this%comm, glb_intrp_find%send_buf(i)%request, ierr)
810 glb_intrp_find%send_buf(i)%flag = .false.
811 end do
812 do i = 1, size(glb_intrp_find%recv_pe)
813 rank = glb_intrp_find%recv_pe(i)
814 call mpi_irecv(this%el_owner0_local(this%n_points_offset_pe_local(rank)+1), &
815 this%n_points_pe_local(rank), &
816 mpi_integer, rank, 0, &
817 this%comm, glb_intrp_find%recv_buf(i)%request, ierr)
818 glb_intrp_find%recv_buf(i)%flag = .false.
819 end do
820 call glb_intrp_find%nbwait_no_op()
821
822 call glb_intrp_find%free()
823
824 !Set up final way of doing communication
825 call send_pe%init()
826 call recv_pe%init()
827 call this%glb_intrp_comm%init_dofs(this%pe_size)
828 do i = 0, (this%pe_size-1)
829 if (this%n_points_pe(i) .gt. 0) then
830 call recv_pe%push(i)
831 point_ids => this%points_at_pe(i)%array()
832 do j = 1, this%n_points_pe(i)
833 call this%glb_intrp_comm%recv_dof(i)%push(point_ids(j))
834 end do
835 end if
836 if (this%n_points_pe_local(i) .gt. 0) then
837 call send_pe%push(i)
838 do j = 1, this%n_points_pe_local(i)
839 call this%glb_intrp_comm%send_dof(i)%push(j+this%n_points_offset_pe_local(i))
840 end do
841 end if
842 end do
843 call this%glb_intrp_comm%init(send_pe, recv_pe,this%comm)
844
845 !Initialize working arrays for evaluation
846 call this%temp_local%init(this%n_points_local)
847 call this%temp%init(this%n_points)
848
849 !Initialize interpolator for local interpolation
850 call this%local_interp%init(this%Xh, this%rst_local, &
851 this%n_points_local)
852
853
854 if (neko_bcknd_device .eq. 1) then
855 call device_memcpy(this%el_owner0, this%el_owner0_d, &
856 this%n_points, host_to_device, sync = .true.)
857 call device_map(this%el_owner0_local, this%el_owner0_local_d, this%n_points_local)
858 call device_memcpy(this%el_owner0_local, this%el_owner0_local_d, &
859 this%n_points_local, host_to_device, sync = .true.)
860 end if
861
862 call this%check_points(this%x%x,this%y%x, this%z%x)
863
864 !Free stuff
865 call send_pe%free()
866 call recv_pe%free()
867 call glb_intrp_find_back%free()
868 call send_pe_find%free()
869 call recv_pe_find%free()
870 call x_t%free()
871 call y_t%free()
872 call z_t%free()
873 call rst_local_cand%free()
874 call resx%free()
875 call resy%free()
876 call resz%free()
877 call res%free()
878 call all_el_candidates%free()
879
880 if (allocated(n_el_cands)) deallocate(n_el_cands)
881 if (allocated(rst_results)) deallocate(rst_results)
882 if (allocated(res_results)) deallocate(res_results)
883 if (allocated(el_owner_results)) deallocate(el_owner_results)
884 call mpi_barrier(this%comm, ierr)
885 time2 = mpi_wtime()
886 write(log_buf, '(A,E15.7)') 'Global interpolation find points done, time (s):', &
887 time2-time_start
888 call neko_log%message(log_buf)
889 call neko_log%end_section()
890 call neko_log%newline()
892
896 subroutine global_interpolation_check_points(this, x, y, z)
897 class(global_interpolation_t), target, intent(inout) :: this
898 real(kind=rp), intent(inout) :: x(:)
899 real(kind=rp), intent(inout) :: y(:)
900 real(kind=rp), intent(inout) :: z(:)
901 integer :: i, j
902 character(len=8000) :: log_buf
903 real(kind=rp) :: xdiff, ydiff, zdiff
904 logical :: isdiff
905 type(vector_t) :: x_check, y_check, z_check
906
907 call x_check%init(this%n_points)
908 call y_check%init(this%n_points)
909 call z_check%init(this%n_points)
910 call this%evaluate(x_check%x, x, on_host=.true.)
911 call this%evaluate(y_check%x, y, on_host=.true.)
912 call this%evaluate(z_check%x, z, on_host=.true.)
913 write(log_buf,'(A)') 'Checking validity of points.'
914 call neko_log%message(log_buf)
915 j = 0
916 do i = 1 , this%n_points
917 ! Check validity of points
918 isdiff = .false.
919 xdiff = x_check%x(i)-this%xyz(1,i)
920 ydiff = y_check%x(i)-this%xyz(2,i)
921 zdiff = z_check%x(i)-this%xyz(3,i)
922 isdiff = norm2(real((/xdiff,ydiff,zdiff/),xp)) > this%tol
923 if (isdiff) then
924 write(*,*) 'Point ', i,'at rank ', this%pe_rank, 'with coordinates: ', &
925 this%xyz(1, i), this%xyz(2, i), this%xyz(3, i), &
926 'Differ from interpolated coords: ', &
927 x_check%x(i), y_check%x(i), z_check%x(i), &
928 'Actual difference: ', &
929 xdiff, ydiff, zdiff, norm2(real((/xdiff,ydiff,zdiff/),xp)),&
930 'Process, element: ', &
931 this%pe_owner(i), this%el_owner0(i)+1, &
932 'Calculated rst: ', &
933 this%rst(1,i), this%rst(2,i), this%rst(3,i)
934 j = j + 1
935 end if
936 end do
937 call x_check%free()
938 call y_check%free()
939 call z_check%free()
940
942
950 subroutine global_interpolation_find_coords(this, x, y, z, n_points)
951 class(global_interpolation_t), intent(inout) :: this
952 integer :: n_points
953 !!Perhaps this should be kind dp
954 !!this is to get around that x,y,z often is 4 dimensional...
955 !!Should maybe add interface for 1d aswell
956 real(kind=rp) :: x(n_points,1,1,1)
957 real(kind=rp) :: y(n_points,1,1,1)
958 real(kind=rp) :: z(n_points,1,1,1)
959 integer :: i
960
961 call this%free_points()
962 call this%free_points_local()
963
964 this%n_points = n_points
965
967
968 !Deepcopy of coordinates
969 do i = 1, n_points
970 this%xyz(1, i) = x(i,1,1,1)
971 this%xyz(2, i) = y(i,1,1,1)
972 this%xyz(3, i) = z(i,1,1,1)
973 end do
974
976
985 subroutine global_interpolation_find_coords1d(this, x, y, z, n_points)
986 class(global_interpolation_t), intent(inout) :: this
987 integer :: n_points
988 real(kind=rp) :: x(n_points)
989 real(kind=rp) :: y(n_points)
990 real(kind=rp) :: z(n_points)
991 integer :: i
992
993 call this%free_points()
994 call this%free_points_local()
995
996 this%n_points = n_points
997
999
1000 !Deepcopy of coordinates
1001 do i = 1, n_points
1002 this%xyz(1, i) = x(i)
1003 this%xyz(2, i) = y(i)
1004 this%xyz(3, i) = z(i)
1005 end do
1006
1008
1010
1011
1013 class(global_interpolation_t) :: this
1014
1015 if (allocated(this%xyz)) deallocate(this%xyz)
1016 if (allocated(this%rst)) deallocate(this%rst)
1017 if (allocated(this%pe_owner)) deallocate(this%pe_owner)
1018 if (allocated(this%el_owner0)) deallocate(this%el_owner0)
1019
1020 allocate(this%pe_owner(this%n_points))
1021 allocate(this%el_owner0(this%n_points))
1022 allocate(this%xyz(3, this%n_points))
1023 allocate(this%rst(3, this%n_points))
1024 if (neko_bcknd_device .eq. 1) then
1025 call device_map(this%el_owner0, this%el_owner0_d, this%n_points)
1026 end if
1027
1029
1036 subroutine global_interpolation_find_xyz(this, xyz, n_points)
1037 class(global_interpolation_t), intent(inout) :: this
1038 integer, intent(in) :: n_points
1039 !!Perhaps this should be kind dp
1040 real(kind=rp), intent(inout) :: xyz(3, n_points)
1041
1042
1043 call this%free_points()
1044 call this%free_points_local()
1045
1046 this%n_points = n_points
1047
1049
1051 call copy(this%xyz, xyz, 3 * n_points)
1052
1054
1055 end subroutine global_interpolation_find_xyz
1056
1065 subroutine global_interpolation_find_and_redist(this, xyz, n_points)
1066 class(global_interpolation_t), intent(inout) :: this
1067 integer, intent(inout) :: n_points
1068 !!Perhaps this should be kind dp
1069 real(kind=rp), allocatable, intent(inout) :: xyz(:,:)
1070
1071
1072 call this%free_points()
1073 call this%free_points_local()
1074
1075
1076 this%n_points = n_points
1078
1080 call copy(this%xyz, xyz, 3 * n_points)
1081
1083 call this%free_points()
1084 this%n_points = this%n_points_local
1085 n_points = this%n_points_local
1087 if (allocated(xyz)) then
1088 deallocate(xyz)
1089 end if
1090 allocate(xyz(3,n_points))
1091
1092 call copy(xyz, this%xyz_local, 3*n_points)
1093 call copy(this%rst, this%rst_local, 3*n_points)
1094 call copy(this%xyz, this%xyz_local, 3*n_points)
1095 this%pe_owner = this%pe_rank
1096 this%el_owner0 = this%el_owner0_local
1097 if (neko_bcknd_device .eq. 1) then
1098 call device_memcpy(this%el_owner0, this%el_owner0_d, &
1099 this%n_points, host_to_device, sync = .true.)
1100 end if
1101 this%all_points_local = .true.
1102
1104
1109 subroutine global_interpolation_evaluate(this, interp_values, field, on_host)
1110 class(global_interpolation_t), target, intent(inout) :: this
1111 real(kind=rp), intent(inout), target :: interp_values(this%n_points)
1112 real(kind=rp), intent(inout), target :: field(this%nelv*this%Xh%lxyz)
1113 logical, intent(in) :: on_host
1114 type(c_ptr) :: interp_d
1115
1116 if (.not. this%all_points_local) then
1117 call this%local_interp%evaluate(this%temp_local%x, this%el_owner0_local, &
1118 field, this%nelv, on_host)
1119 if (neko_bcknd_device .eq. 1 .and. .not. on_host) then
1120 call device_memcpy(this%temp_local%x, this%temp_local%x_d, &
1121 this%n_points_local, device_to_host, .true.)
1122 end if
1123 interp_values = 0.0_rp
1124 call this%glb_intrp_comm%sendrecv(this%temp_local%x, interp_values, &
1125 this%n_points_local, this%n_points)
1126 if (neko_bcknd_device .eq. 1 .and. .not. on_host) then
1127 interp_d = device_get_ptr(interp_values)
1128 call device_memcpy(interp_values, interp_d, &
1129 this%n_points, host_to_device, .false.)
1130 end if
1131 else
1132 call this%local_interp%evaluate(interp_values, this%el_owner0_local, &
1133 field, this%nelv, on_host)
1134 end if
1135
1136 end subroutine global_interpolation_evaluate
1137
1138
1149 function rst_cmp(rst1, rst2,res1, res2, tol) result(rst2_better)
1150 real(kind=rp) :: rst1(3), res1(3)
1151 real(kind=rp) :: rst2(3), res2(3)
1152 real(kind=rp) :: tol
1153 logical :: rst2_better
1154 !If rst1 is invalid and rst2 is valid, take rst2
1155 ! If both invalidl, take smallest residual
1156 if (abs(rst1(1)) .gt. 1.0_xp+tol .or. &
1157 abs(rst1(2)) .gt. 1.0_xp+tol .or. &
1158 abs(rst1(3)) .gt. 1.0_xp+tol) then
1159 if (abs(rst2(1)) .le. 1.0_xp+tol .and. &
1160 abs(rst2(2)) .le. 1.0_xp+tol .and. &
1161 abs(rst2(3)) .le. 1.0_xp+tol) then
1162 rst2_better = .true.
1163 else
1164 rst2_better = (norm2(real(res2,xp)) .lt. norm2(real(res1,xp)))
1165 end if
1166 else
1168 rst2_better = (norm2(real(res2,xp)) .lt. norm2(real(res1,xp)) .and.&
1169 abs(rst2(1)) .le. 1.0_xp+tol .and. &
1170 abs(rst2(2)) .le. 1.0_xp+tol .and. &
1171 abs(rst2(3)) .le. 1.0_xp+tol)
1172 end if
1173 end function rst_cmp
1174
1175end module global_interpolation
double real
Deassociate a Fortran array from a device pointer.
Definition device.F90:95
Return the device pointer for an associated Fortran array.
Definition device.F90:101
Map a Fortran array to a device (allocate and associate)
Definition device.F90:77
Copy data between host and device (or device and device)
Definition device.F90:71
Implements aabb_pe_finder given a dofmap.
integer, parameter, public glob_map_size
Minimum number of total boxes in the aabb tree.
Implements cartesian_pe_finder given a dofmap.
Definition comm.F90:1
integer, public pe_size
MPI size of communicator.
Definition comm.F90:58
integer, public pe_rank
MPI rank.
Definition comm.F90:55
type(mpi_comm), public neko_comm
MPI communicator.
Definition comm.F90:42
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:219
integer, parameter, public device_to_host
Definition device.F90:47
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
Defines a field.
Definition field.f90:34
Defines global interpolation communication Based on the MPI based gather-scatter kernel.
Implements global_interpolation given a dofmap.
subroutine global_interpolation_free(this)
Destructor.
subroutine global_interpolation_find_coords1d(this, x, y, z, n_points)
Finds the corresponding r,s,t coordinates in the correct global element as well as which process that...
subroutine global_interpolation_free_points(this)
Destructor for point arrays.
subroutine global_interpolation_init_point_arrays(this)
subroutine global_interpolation_evaluate(this, interp_values, field, on_host)
Evalute the interpolated value in the points given a field.
subroutine global_interpolation_free_points_local(this)
subroutine global_interpolation_find_coords(this, x, y, z, n_points)
Finds the corresponding r,s,t coordinates in the correct global element as well as which process that...
logical function rst_cmp(rst1, rst2, res1, res2, tol)
Compares two sets of rst coordinates and checks whether rst2 is better than rst1 given a tolerance re...
subroutine global_interpolation_init_dof(this, dof, comm, tol, pad)
Initialize the global interpolation object on a dofmap.
subroutine global_interpolation_init_xyz(this, x, y, z, gdim, nelv, xh, comm, tol, pad)
Initialize the global interpolation object on a set of coordinates.
subroutine global_interpolation_find_xyz(this, xyz, n_points)
Finds the corresponding r,s,t coordinates in the correct global element as well as which process that...
subroutine global_interpolation_find_and_redist(this, xyz, n_points)
Finds the corresponding r,s,t coordinates and redistributes the points to the owning rank in the corr...
subroutine global_interpolation_find_common(this)
Common routine for finding the points.
subroutine global_interpolation_check_points(this, x, y, z)
Check the points for validity This is used to check that the points are valid and that the interpolat...
Routines to obtain interpolated values on a set of points with known rst coordinates in elements loca...
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
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:249
real(kind=rp), parameter, public neko_eps
Machine epsilon .
Definition math.f90:69
Defines a matrix.
Definition matrix.f90:34
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public xp
Definition num_types.f90:14
integer, parameter, public dp
Definition num_types.f90:9
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Implements pe_finder given a dofmap.
Definition pe_finder.f90:35
Defines a function space.
Definition space.f90:34
Implements a dynamic stack ADT.
Definition stack.f90:35
Defines structs that are used... Dont know if we should keep it though.
Definition structs.f90:2
Utilities.
Definition utils.f90:35
Defines a vector.
Definition vector.f90:34
Implements global interpolation for arbitrary points in the domain.
Implements global interpolation for arbitrary points in the domain.
Implements global interpolation for arbitrary points in the domain.
Base type for element finder providing element candidates for a given point in the domain.
Definition el_finder.f90:43
Global interpolation communication method.
Implements global interpolation for arbitrary points in the domain.
Type to compute local element (rst) coordinates for a gives set points in physical (xyz) space on a S...
Interpolation on a set of points with known rst coordinates in elements local to this process....
Implements global interpolation for arbitrary points in the domain.
Definition pe_finder.f90:45
The function space for the SEM solution fields.
Definition space.f90:63
Integer based stack.
Definition stack.f90:63
Pointer to array.
Definition structs.f90:14
#define max(a, b)
Definition tensor.cu:40