Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
bc.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!
34module bc
36 use num_types, only : rp
37 use device, only : host_to_device, device_memcpy, &
39 use iso_c_binding, only: c_associated
40 use dofmap, only : dofmap_t
41 use coefs, only : coef_t
42 use space, only : space_t
44 use facet_zone, only : facet_zone_t
45 use stack, only : stack_i4t2_t
46 use tuple, only : tuple_i4_t
47 use field, only : field_t
48 use gs_ops, only : gs_op_add
49 use math, only : relcmp
51 use, intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr
52 use json_module, only : json_file
53 use time_state, only : time_state_t
54 use field, only : field_t
55 use file, only : file_t
56
57 implicit none
58 private
59
61 type, public, abstract :: bc_t
63 integer, allocatable :: msk(:)
65 integer, allocatable :: facet(:)
67 type(dofmap_t), pointer :: dof => null()
69 type(coef_t), pointer :: coef => null()
71 type(mesh_t), pointer :: msh => null()
73 type(space_t), pointer :: xh => null()
75 type(stack_i4t2_t) :: marked_facet
77 type(c_ptr) :: msk_d = c_null_ptr
79 type(c_ptr) :: facet_d = c_null_ptr
84 logical :: strong = .true.
87 logical :: updated = .false.
88 contains
90 procedure, pass(this) :: init_base => bc_init_base
92 procedure, pass(this) :: free_base => bc_free_base
94 procedure, pass(this) :: mark_facet => bc_mark_facet
96 procedure, pass(this) :: mark_facets => bc_mark_facets
98 procedure, pass(this) :: mark_zone => bc_mark_zone
101 procedure, pass(this) :: finalize_base => bc_finalize_base
102
105 procedure, pass(this) :: apply_scalar_generic => bc_apply_scalar_generic
108 procedure, pass(this) :: apply_vector_generic => bc_apply_vector_generic
110 procedure, pass(this) :: debug_mask_ => bc_debug_mask
112 procedure(bc_apply_scalar), pass(this), deferred :: apply_scalar
114 procedure(bc_apply_vector), pass(this), deferred :: apply_vector
116 procedure(bc_apply_scalar_dev), pass(this), deferred :: apply_scalar_dev
118 procedure(bc_apply_vector_dev), pass(this), deferred :: apply_vector_dev
120 procedure(bc_destructor), pass(this), deferred :: free
122 procedure(bc_constructor), pass(this), deferred :: init
124 procedure(bc_finalize), pass(this), deferred :: finalize
125 end type bc_t
126
128 type, public :: bc_ptr_t
129 class(bc_t), pointer :: ptr => null()
130 end type bc_ptr_t
131
132 ! Helper type to have an array of polymorphic bc_t objects.
133 type, public :: bc_alloc_t
134 class(bc_t), allocatable :: obj
135 end type bc_alloc_t
136
137
138 abstract interface
139
140 subroutine bc_constructor(this, coef, json)
141 import :: bc_t, coef_t, json_file
142 class(bc_t), intent(inout), target :: this
143 type(coef_t), target, intent(in) :: coef
144 type(json_file), intent(inout) :: json
145 end subroutine bc_constructor
146 end interface
147
148 abstract interface
149
150 subroutine bc_destructor(this)
151 import :: bc_t
152 class(bc_t), intent(inout), target :: this
153 end subroutine bc_destructor
154 end interface
155
156 abstract interface
157
158 subroutine bc_finalize(this, only_facets)
159 import :: bc_t
160 class(bc_t), intent(inout), target :: this
161 logical, optional, intent(in) :: only_facets
162 end subroutine bc_finalize
163 end interface
164
165 abstract interface
166
171 subroutine bc_apply_scalar(this, x, n, time, strong)
172 import :: bc_t, time_state_t
173 import :: rp
174 class(bc_t), intent(inout) :: this
175 integer, intent(in) :: n
176 real(kind=rp), intent(inout), dimension(n) :: x
177 type(time_state_t), intent(in), optional :: time
178 logical, intent(in), optional :: strong
179 end subroutine bc_apply_scalar
180 end interface
181
182 abstract interface
183
191 subroutine bc_apply_vector(this, x, y, z, n, time, strong)
192 import :: bc_t, time_state_t
193 import :: rp
194 class(bc_t), intent(inout) :: this
195 integer, intent(in) :: n
196 real(kind=rp), intent(inout), dimension(n) :: x
197 real(kind=rp), intent(inout), dimension(n) :: y
198 real(kind=rp), intent(inout), dimension(n) :: z
199 type(time_state_t), intent(in), optional :: time
200 logical, intent(in), optional :: strong
201 end subroutine bc_apply_vector
202 end interface
203
204 abstract interface
205
210 subroutine bc_apply_scalar_dev(this, x_d, time, strong, strm)
211 import :: c_ptr
212 import :: bc_t, time_state_t
213 import :: rp
214 class(bc_t), intent(inout), target :: this
215 type(c_ptr), intent(inout) :: x_d
216 type(time_state_t), intent(in), optional :: time
217 logical, intent(in), optional :: strong
218 type(c_ptr), intent(inout) :: strm
219 end subroutine bc_apply_scalar_dev
220 end interface
221
222 abstract interface
223
230 subroutine bc_apply_vector_dev(this, x_d, y_d, z_d, time, strong, strm)
231 import :: c_ptr, bc_t, time_state_t
232 import :: rp
233 class(bc_t), intent(inout), target :: this
234 type(c_ptr), intent(inout) :: x_d
235 type(c_ptr), intent(inout) :: y_d
236 type(c_ptr), intent(inout) :: z_d
237 type(time_state_t), intent(in), optional :: time
238 logical, intent(in), optional :: strong
239 type(c_ptr), intent(inout) :: strm
240 end subroutine bc_apply_vector_dev
241 end interface
242
243contains
244
247 subroutine bc_init_base(this, coef)
248 class(bc_t), intent(inout) :: this
249 type(coef_t), target, intent(in) :: coef
250
251 call this%free_base
252
253 this%dof => coef%dof
254 this%coef => coef
255 this%Xh => this%dof%Xh
256 this%msh => this%dof%msh
257
258 call this%marked_facet%init()
259
260 end subroutine bc_init_base
261
263 subroutine bc_free_base(this)
264 class(bc_t), intent(inout) :: this
265
266 call this%marked_facet%free()
267
268 nullify(this%Xh)
269 nullify(this%msh)
270 nullify(this%dof)
271 nullify(this%coef)
272
273 if (allocated(this%msk)) then
274 deallocate(this%msk)
275 end if
276
277 if (allocated(this%facet)) then
278 deallocate(this%facet)
279 end if
280
281 if (c_associated(this%msk_d)) then
282 call device_free(this%msk_d)
283 this%msk_d = c_null_ptr
284 end if
285
286 if (c_associated(this%facet_d)) then
287 call device_free(this%facet_d)
288 this%facet_d = c_null_ptr
289 end if
290
291 end subroutine bc_free_base
292
301 subroutine bc_apply_vector_generic(this, x, y, z, time, strong, strm)
302 class(bc_t), intent(inout) :: this
303 type(field_t), intent(inout) :: x
304 type(field_t), intent(inout) :: y
305 type(field_t), intent(inout) :: z
306 type(time_state_t), intent(in), optional :: time
307 logical, intent(in), optional :: strong
308 type(c_ptr), intent(inout), optional :: strm
309 type(c_ptr) :: strm_
310 integer :: n
311 character(len=256) :: msg
312
313 ! Get the size of the fields
314 n = x%size()
315
316 ! Ensure all fields are the same size
317 if (y%size() .ne. n .or. z%size() .ne. n) then
318 msg = "Fields x, y, z must have the same size in " // &
319 "bc_list_apply_vector_field"
320 call neko_error(trim(msg))
321 end if
322
323 if (neko_bcknd_device .eq. 1) then
324
325 if (present(strm)) then
326 strm_ = strm
327 else
328 strm_ = glb_cmd_queue
329 end if
330
331 call this%apply_vector_dev(x%x_d, y%x_d, z%x_d, time = time, &
332 strong = strong, strm = strm_)
333 else
334 call this%apply_vector(x%x, y%x, z%x, n, time = time, strong = strong)
335 end if
336
337 end subroutine bc_apply_vector_generic
338
345 subroutine bc_apply_scalar_generic(this, x, time, strong, strm)
346 class(bc_t), intent(inout) :: this
347 type(field_t), intent(inout) :: x
348 type(time_state_t), intent(in), optional :: time
349 logical, intent(in), optional :: strong
350 type(c_ptr), intent(inout), optional :: strm
351 type(c_ptr) :: strm_
352 integer :: n
353
354 ! Get the size of the field
355 n = x%size()
356
357 if (neko_bcknd_device .eq. 1) then
358
359 if (present(strm)) then
360 strm_ = strm
361 else
362 strm_ = glb_cmd_queue
363 end if
364
365 call this%apply_scalar_dev(x%x_d, time = time, strong = strong, &
366 strm = strm_)
367 else
368 call this%apply_scalar(x%x, n, time = time)
369 end if
370
371 end subroutine bc_apply_scalar_generic
372
376 subroutine bc_mark_facet(this, facet, el)
377 class(bc_t), intent(inout) :: this
378 integer, intent(in) :: facet
379 integer, intent(in) :: el
380 type(tuple_i4_t) :: t
381
382 t%x = [facet, el]
383 call this%marked_facet%push(t)
384
385 end subroutine bc_mark_facet
386
389 subroutine bc_mark_facets(this, facet_list)
390 class(bc_t), intent(inout) :: this
391 type(stack_i4t2_t), intent(inout) :: facet_list
392 type(tuple_i4_t), pointer :: fp(:)
393 integer :: i
394
395 fp => facet_list%array()
396 do i = 1, facet_list%size()
397 call this%marked_facet%push(fp(i))
398 end do
399
400 end subroutine bc_mark_facets
401
404 subroutine bc_mark_zone(this, bc_zone)
405 class(bc_t), intent(inout) :: this
406 class(facet_zone_t), intent(in) :: bc_zone
407 integer :: i
408 do i = 1, bc_zone%size
409 call this%marked_facet%push(bc_zone%facet_el(i))
410 end do
411 end subroutine bc_mark_zone
412
420 subroutine bc_finalize_base(this, only_facets)
421 class(bc_t), target, intent(inout) :: this
422 logical, optional, intent(in) :: only_facets
423 type(tuple_i4_t), pointer :: bfp(:)
424 type(tuple_i4_t) :: bc_facet
425 type(field_t) :: test_field
426 integer :: facet_size, facet, el
427 logical :: only_facet = .false.
428 integer :: i, j, k, l, msk_c
429 integer :: lx, ly, lz, n
430 lx = this%Xh%lx
431 ly = this%Xh%ly
432 lz = this%Xh%lz
433 if ( present(only_facets)) then
434 only_facet = only_facets
435 else
436 only_facet = .false.
437 end if
439
440 ! Note we assume that lx = ly = lz
441 facet_size = lx**2
442 allocate(this%msk(0:facet_size * this%marked_facet%size()))
443 allocate(this%facet(0:facet_size * this%marked_facet%size()))
444
445 msk_c = 0
446 bfp => this%marked_facet%array()
447
448 ! Loop through each (facet, element) id tuple
449 ! Then loop over all the nodes of the face and compute their linear index
450 ! This index goes into this%msk, whereas the corresponding face id goes into
451 ! this%facet
452 do i = 1, this%marked_facet%size()
453 bc_facet = bfp(i)
454 facet = bc_facet%x(1)
455 el = bc_facet%x(2)
456 select case (facet)
457 case (1)
458 do l = 1, lz
459 do k = 1, ly
460 msk_c = msk_c + 1
461 this%msk(msk_c) = linear_index(1, k, l, el, lx, ly, lz)
462 this%facet(msk_c) = 1
463 end do
464 end do
465 case (2)
466 do l = 1, lz
467 do k = 1, ly
468 msk_c = msk_c + 1
469 this%msk(msk_c) = linear_index(lx, k, l, el, lx, ly, lz)
470 this%facet(msk_c) = 2
471 end do
472 end do
473 case (3)
474 do l = 1, lz
475 do j = 1, lx
476 msk_c = msk_c + 1
477 this%msk(msk_c) = linear_index(j, 1, l, el, lx, ly, lz)
478 this%facet(msk_c) = 3
479 end do
480 end do
481 case (4)
482 do l = 1, lz
483 do j = 1, lx
484 msk_c = msk_c + 1
485 this%msk(msk_c) = linear_index(j, ly, l, el, lx, ly, lz)
486 this%facet(msk_c) = 4
487 end do
488 end do
489 case (5)
490 do k = 1, ly
491 do j = 1, lx
492 msk_c = msk_c + 1
493 this%msk(msk_c) = linear_index(j, k, 1, el, lx, ly, lz)
494 this%facet(msk_c) = 5
495 end do
496 end do
497 case (6)
498 do k = 1, ly
499 do j = 1, lx
500 msk_c = msk_c + 1
501 this%msk(msk_c) = linear_index(j, k, lz, el, lx, ly, lz)
502 this%facet(msk_c) = 6
503 end do
504 end do
505 end select
506 end do
507 this%facet(0) = msk_c
508 if (neko_bcknd_device .eq. 1) then
509 !Observe the facet_mask is junk if only_facet is false
510 n = msk_c + 1
511 call device_map(this%facet, this%facet_d, n)
512 call device_memcpy(this%facet, this%facet_d, n, &
513 host_to_device, sync = .true.)
514 end if
515 if ( .not. only_facet) then
516 !Makes check for points not on facet that should have bc applied
517 call test_field%init(this%dof)
518
519 n = test_field%size()
520 test_field%x = 0.0_rp
521 !Apply this bc once
522 do i = 1, msk_c
523 test_field%x(this%msk(i),1,1,1) = 1.0
524 end do
525 if (neko_bcknd_device .eq. 1) then
526 call device_memcpy(test_field%x, test_field%x_d, n, &
527 host_to_device, sync=.true.)
528 end if
529 !Check if some point that was not zeroed was zeroed on another element
530 call this%coef%gs_h%op(test_field, gs_op_add)
531 if (neko_bcknd_device .eq. 1) then
532 call device_memcpy(test_field%x, test_field%x_d, n, &
533 device_to_host, sync=.true.)
534 end if
535 msk_c = 0
536 do i = 1, this%dof%size()
537 if (test_field%x(i,1,1,1) .gt. 0.5) then
538 msk_c = msk_c + 1
539 end if
540 end do
541 !Allocate new mask
542 deallocate(this%msk)
543 allocate(this%msk(0:msk_c))
544 j = 1
545 do i = 1, this%dof%size()
546 if (test_field%x(i,1,1,1) .gt. 0.5) then
547 this%msk(j) = i
548 j = j + 1
549 end if
550 end do
551
552 call test_field%free()
553 end if
554
555 this%msk(0) = msk_c
556
557 if (neko_bcknd_device .eq. 1) then
558 n = msk_c + 1
559 call device_map(this%msk, this%msk_d, n)
560 call device_memcpy(this%msk, this%msk_d, n, &
561 host_to_device, sync = .true.)
562 end if
563
564 end subroutine bc_finalize_base
565
569 subroutine bc_debug_mask(this, file_name)
570 class(bc_t), intent(inout) :: this
571 character(len=*), intent(in) :: file_name
572 type(field_t) :: bdry_field
573 integer:: i, m, k
574 type(file_t) :: dump_file
575
576 call bdry_field%init(this%coef%dof, 'bdry')
577 m = this%msk(0)
578 do i = 1, m
579 k = this%msk(i)
580 bdry_field%x(k,1,1,1) = 1.0_rp
581 end do
582 call dump_file%init(file_name)
583 call dump_file%write(bdry_field)
584
585 end subroutine bc_debug_mask
586end module bc
Apply the boundary condition to a scalar field on the device.
Definition bc.f90:210
Apply the boundary condition to a scalar field.
Definition bc.f90:171
Apply the boundary condition to a vector field on the device.
Definition bc.f90:230
Apply the boundary condition to a vector field.
Definition bc.f90:191
Constructor.
Definition bc.f90:140
Destructor.
Definition bc.f90:150
Finalize by building the mask and facet arrays.
Definition bc.f90:158
Map a Fortran array to a device (allocate and associate)
Definition device.F90:72
Copy data between host and device (or device and device)
Definition device.F90:66
Defines a boundary condition.
Definition bc.f90:34
subroutine bc_mark_zone(this, bc_zone)
Mark all facets from a zone.
Definition bc.f90:405
subroutine bc_free_base(this)
Destructor for the base type, bc_t.
Definition bc.f90:264
subroutine bc_init_base(this, coef)
Constructor.
Definition bc.f90:248
subroutine bc_apply_scalar_generic(this, x, time, strong, strm)
Apply the boundary condition to a scalar field. Dispatches to the CPU or the device version.
Definition bc.f90:346
subroutine bc_finalize_base(this, only_facets)
Finalize the construction of the bc by populting the msk and facet arrays.
Definition bc.f90:421
subroutine bc_apply_vector_generic(this, x, y, z, time, strong, strm)
Apply the boundary condition to a vector field. Dispatches to the CPU or the device version.
Definition bc.f90:302
subroutine bc_debug_mask(this, file_name)
Write a field showing the mask of the bc.
Definition bc.f90:570
subroutine bc_mark_facet(this, facet, el)
Mark facet on element el as part of the boundary condition.
Definition bc.f90:377
subroutine bc_mark_facets(this, facet_list)
Mark all facets from a (facet, el) tuple list.
Definition bc.f90:390
Coefficients.
Definition coef.f90:34
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:214
integer, parameter, public device_to_host
Definition device.F90:47
type(c_ptr), bind(C), public glb_cmd_queue
Global command queue.
Definition device.F90:51
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
Defines a zone as a subset of facets in a mesh.
Defines a field.
Definition field.f90:34
Module for file I/O operations.
Definition file.f90:34
Defines Gather-scatter operations.
Definition gs_ops.f90:34
integer, parameter, public gs_op_add
Definition gs_ops.f90:36
Definition math.f90:60
Defines a mesh.
Definition mesh.f90:34
integer, parameter, public neko_msh_max_zlbls
Max num. zone labels.
Definition mesh.f90:62
integer, parameter, public neko_msh_max_zlbl_len
Max length of a zone label.
Definition mesh.f90:64
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Defines a function space.
Definition space.f90:34
Implements a dynamic stack ADT.
Definition stack.f90:35
Module with things related to the simulation time.
Implements a n-tuple.
Definition tuple.f90:34
Utilities.
Definition utils.f90:35
character(len=100) function, dimension(:), allocatable, public split_string(string, delimiter)
Split a string based on delimiter (tokenizer) OBS: very hacky, this should really be improved,...
Definition utils.f90:137
pure integer function, public linear_index(i, j, k, l, lx, ly, lz)
Compute the address of a (i,j,k,l) array with sizes (1:lx, 1:ly, 1:lz, :)
Definition utils.f90:175
Pointer to a `bc_t`.
Definition bc.f90:128
Base type for a boundary condition.
Definition bc.f90:61
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
A wrapper around a polymorphic generic_file_t that handles its init. This is essentially a factory fo...
Definition file.f90:55
The function space for the SEM solution fields.
Definition space.f90:62
Integer 2-tuple based stack.
Definition stack.f90:84
A struct that contains all info about the time, expand as needed.
Integer based 2-tuple.
Definition tuple.f90:51