Neko 0.9.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-2024, 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
35 use neko_config
36 use num_types
37 use device
38 use dofmap, only : dofmap_t
39 use coefs, only : coef_t
40 use space, only : space_t
42 use facet_zone, only : facet_zone_t
43 use stack, only : stack_i4t2_t
44 use tuple, only : tuple_i4_t
46 use, intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr
47 implicit none
48 private
49
51 type, public, abstract :: bc_t
53 integer, allocatable :: msk(:)
55 integer, allocatable :: facet(:)
57 type(dofmap_t), pointer :: dof
59 type(coef_t), pointer :: coef
61 type(mesh_t), pointer :: msh
63 type(space_t), pointer :: xh
65 type(stack_i4t2_t) :: marked_facet
67 type(c_ptr) :: msk_d = c_null_ptr
69 type(c_ptr) :: facet_d = c_null_ptr
70 contains
72 procedure, pass(this) :: init_base => bc_init_base
74 procedure, pass(this) :: free_base => bc_free_base
76 procedure, pass(this) :: mark_facet => bc_mark_facet
78 procedure, pass(this) :: mark_facets => bc_mark_facets
80 procedure, pass(this) :: mark_zones_from_list => bc_mark_zones_from_list
82 procedure, pass(this) :: mark_zone => bc_mark_zone
85 procedure, pass(this) :: finalize => bc_finalize
87 procedure(bc_apply_scalar), pass(this), deferred :: apply_scalar
89 procedure(bc_apply_vector), pass(this), deferred :: apply_vector
91 procedure(bc_apply_scalar_dev), pass(this), deferred :: apply_scalar_dev
93 procedure(bc_apply_vector_dev), pass(this), deferred :: apply_vector_dev
95 procedure(bc_destructor), pass(this), deferred :: free
96 end type bc_t
97
99 type, private :: bcp_t
100 class(bc_t), pointer :: bcp
101 end type bcp_t
102
104 type, public :: bc_list_t
105 type(bcp_t), allocatable :: bc(:)
107 integer :: n
109 integer :: size
110 end type bc_list_t
111
112 abstract interface
113
118 subroutine bc_apply_scalar(this, x, n, t, tstep)
119 import :: bc_t
120 import :: rp
121 class(bc_t), intent(inout) :: this
122 integer, intent(in) :: n
123 real(kind=rp), intent(inout), dimension(n) :: x
124 real(kind=rp), intent(in), optional :: t
125 integer, intent(in), optional :: tstep
126 end subroutine bc_apply_scalar
127 end interface
128
129 abstract interface
130
137 subroutine bc_apply_vector(this, x, y, z, n, t, tstep)
138 import :: bc_t
139 import :: rp
140 class(bc_t), intent(inout) :: this
141 integer, intent(in) :: n
142 real(kind=rp), intent(inout), dimension(n) :: x
143 real(kind=rp), intent(inout), dimension(n) :: y
144 real(kind=rp), intent(inout), dimension(n) :: z
145 real(kind=rp), intent(in), optional :: t
146 integer, intent(in), optional :: tstep
147 end subroutine bc_apply_vector
148 end interface
149
150 abstract interface
151
152 subroutine bc_destructor(this)
153 import :: bc_t
154 class(bc_t), intent(inout), target :: this
155 end subroutine bc_destructor
156 end interface
157
158 abstract interface
159
161 subroutine bc_apply_scalar_dev(this, x_d, t, tstep)
162 import :: c_ptr
163 import :: bc_t
164 import :: rp
165 class(bc_t), intent(inout), target :: this
166 type(c_ptr) :: x_d
167 real(kind=rp), intent(in), optional :: t
168 integer, intent(in), optional :: tstep
169 end subroutine bc_apply_scalar_dev
170 end interface
171
172 abstract interface
173
177 subroutine bc_apply_vector_dev(this, x_d, y_d, z_d, t, tstep)
178 import :: c_ptr
179 import :: bc_t
180 import :: rp
181 class(bc_t), intent(inout), target :: this
182 type(c_ptr) :: x_d
183 type(c_ptr) :: y_d
184 type(c_ptr) :: z_d
185 real(kind=rp), intent(in), optional :: t
186 integer, intent(in), optional :: tstep
187 end subroutine bc_apply_vector_dev
188 end interface
189
192 end interface bc_list_apply
193
196
197contains
198
201 subroutine bc_init_base(this, coef)
202 class(bc_t), intent(inout) :: this
203 type(coef_t), target, intent(in) :: coef
204
205 call this%free_base
206
207 this%dof => coef%dof
208 this%coef => coef
209 this%Xh => this%dof%Xh
210 this%msh => this%dof%msh
211
212 call this%marked_facet%init()
213
214 end subroutine bc_init_base
215
217 subroutine bc_free_base(this)
218 class(bc_t), intent(inout) :: this
219
220 call this%marked_facet%free()
221
222 nullify(this%Xh)
223 nullify(this%msh)
224 nullify(this%dof)
225 nullify(this%coef)
226
227 if (allocated(this%msk)) then
228 deallocate(this%msk)
229 end if
230
231 if (allocated(this%facet)) then
232 deallocate(this%facet)
233 end if
234
235 if (c_associated(this%msk_d)) then
236 call device_free(this%msk_d)
237 this%msk_d = c_null_ptr
238 end if
239
240 if (c_associated(this%facet_d)) then
241 call device_free(this%facet_d)
242 this%facet_d = c_null_ptr
243 end if
244
245 end subroutine bc_free_base
246
250 subroutine bc_mark_facet(this, facet, el)
251 class(bc_t), intent(inout) :: this
252 integer, intent(in) :: facet
253 integer, intent(in) :: el
254 type(tuple_i4_t) :: t
255
256 t%x = (/facet, el/)
257 call this%marked_facet%push(t)
258
259 end subroutine bc_mark_facet
260
263 subroutine bc_mark_facets(this, facet_list)
264 class(bc_t), intent(inout) :: this
265 type(stack_i4t2_t), intent(inout) :: facet_list
266 type(tuple_i4_t), pointer :: fp(:)
267 integer :: i
268
269 fp => facet_list%array()
270 do i = 1, facet_list%size()
271 call this%marked_facet%push(fp(i))
272 end do
273
274 end subroutine bc_mark_facets
275
278 subroutine bc_mark_zone(this, bc_zone)
279 class(bc_t), intent(inout) :: this
280 class(facet_zone_t), intent(inout) :: bc_zone
281 integer :: i
282 do i = 1, bc_zone%size
283 call this%marked_facet%push(bc_zone%facet_el(i))
284 end do
285 end subroutine bc_mark_zone
286
293 subroutine bc_mark_zones_from_list(this, bc_zones, bc_key, bc_labels)
294 class(bc_t), intent(inout) :: this
295 class(facet_zone_t), intent(inout) :: bc_zones(:)
296 character(len=*) :: bc_key
297 character(len=100), allocatable :: split_key(:)
298 character(len=NEKO_MSH_MAX_ZLBL_LEN) :: bc_labels(:)
299 integer :: i, j, k, l, msh_bc_type
300
301 msh_bc_type = 0
302 if(trim(bc_key) .eq. 'o' .or. trim(bc_key) .eq. 'on' &
303 .or. trim(bc_key) .eq. 'o+dong' .or. trim(bc_key) .eq. 'on+dong') then
304 msh_bc_type = 1
305 else if(trim(bc_key) .eq. 'd_pres') then
306 msh_bc_type = 1
307 else if(trim(bc_key) .eq. 'w') then
308 msh_bc_type = 2
309 else if(trim(bc_key) .eq. 'v') then
310 msh_bc_type = 2
311 else if(trim(bc_key) .eq. 'd_vel_u') then
312 msh_bc_type = 2
313 else if(trim(bc_key) .eq. 'd_vel_v') then
314 msh_bc_type = 2
315 else if(trim(bc_key) .eq. 'd_vel_w') then
316 msh_bc_type = 2
317 else if(trim(bc_key) .eq. 'sym') then
318 msh_bc_type = 2
319 else if(trim(bc_key) .eq. 'sh') then
320 msh_bc_type = 2
321 else if(trim(bc_key) .eq. 'wm') then
322 msh_bc_type = 2
323 end if
324
325 do i = 1, size(bc_labels)
326 !Check if several bcs are defined for this zone
327 !bcs are seperated by /, but we could use something else
328 if (index(trim(bc_labels(i)), '/') .eq. 0) then
329 if (trim(bc_key) .eq. trim(bc_labels(i))) then
330 call bc_mark_zone(this, bc_zones(i))
331 ! Loop across all faces in the mesh
332 do j = 1,this%msh%nelv
333 do k = 1, 2 * this%msh%gdim
334 if (this%msh%facet_type(k,j) .eq. -i) then
335 this%msh%facet_type(k,j) = msh_bc_type
336 end if
337 end do
338 end do
339 end if
340 else
341 split_key = split_string(trim(bc_labels(i)),'/')
342 do l = 1, size(split_key)
343 if (trim(split_key(l)) .eq. trim(bc_key)) then
344 call bc_mark_zone(this, bc_zones(i))
345 ! Loop across all faces in the mesh
346 do j = 1,this%msh%nelv
347 do k = 1, 2 * this%msh%gdim
348 if (this%msh%facet_type(k,j) .eq. -i) then
349 this%msh%facet_type(k,j) = msh_bc_type
350 end if
351 end do
352 end do
353 end if
354 end do
355 end if
356 end do
357 end subroutine bc_mark_zones_from_list
358
359
363 subroutine bc_finalize(this)
364 class(bc_t), target, intent(inout) :: this
365 type(tuple_i4_t), pointer :: bfp(:)
366 type(tuple_i4_t) :: bc_facet
367 integer :: facet_size, facet, el
368 integer :: i, j, k, l, msk_c
369 integer :: lx, ly, lz, n
370
371 lx = this%Xh%lx
372 ly = this%Xh%ly
373 lz = this%Xh%lz
374
376
377 ! Note we assume that lx = ly = lz
378 facet_size = lx**2
379 allocate(this%msk(0:facet_size * this%marked_facet%size()))
380 allocate(this%facet(0:facet_size * this%marked_facet%size()))
381
382 msk_c = 0
383 bfp => this%marked_facet%array()
384
385 ! Loop through each (facet, element) id tuple
386 ! Then loop over all the nodes of the face and compute their linear index
387 ! This index goes into this%msk, whereas the corresponding face id goes into
388 ! this%facet
389 do i = 1, this%marked_facet%size()
390 bc_facet = bfp(i)
391 facet = bc_facet%x(1)
392 el = bc_facet%x(2)
393 select case (facet)
394 case (1)
395 do l = 1, lz
396 do k = 1, ly
397 msk_c = msk_c + 1
398 this%msk(msk_c) = linear_index(1,k,l,el,lx,ly,lz)
399 this%facet(msk_c) = 1
400 end do
401 end do
402 case (2)
403 do l = 1, lz
404 do k = 1, ly
405 msk_c = msk_c + 1
406 this%msk(msk_c) = linear_index(lx,k,l,el,lx,ly,lz)
407 this%facet(msk_c) = 2
408 end do
409 end do
410 case(3)
411 do l = 1, lz
412 do j = 1, lx
413 msk_c = msk_c + 1
414 this%msk(msk_c) = linear_index(j,1,l,el,lx,ly,lz)
415 this%facet(msk_c) = 3
416 end do
417 end do
418 case(4)
419 do l = 1, lz
420 do j = 1, lx
421 msk_c = msk_c + 1
422 this%msk(msk_c) = linear_index(j,ly,l,el,lx,ly,lz)
423 this%facet(msk_c) = 4
424 end do
425 end do
426 case(5)
427 do k = 1, ly
428 do j = 1, lx
429 msk_c = msk_c + 1
430 this%msk(msk_c) = linear_index(j,k,1,el,lx,ly,lz)
431 this%facet(msk_c) = 5
432 end do
433 end do
434 case(6)
435 do k = 1, ly
436 do j = 1, lx
437 msk_c = msk_c + 1
438 this%msk(msk_c) = linear_index(j,k,lz,el,lx,ly,lz)
439 this%facet(msk_c) = 6
440 end do
441 end do
442 end select
443 end do
444
445 this%msk(0) = msk_c
446 this%facet(0) = msk_c
447
448 if (neko_bcknd_device .eq. 1) then
449 n = facet_size * this%marked_facet%size() + 1
450 call device_map(this%msk, this%msk_d, n)
451 call device_map(this%facet, this%facet_d, n)
452
453 call device_memcpy(this%msk, this%msk_d, n, &
454 host_to_device, sync=.false.)
455 call device_memcpy(this%facet, this%facet_d, n, &
456 host_to_device, sync=.false.)
457 end if
458
459 end subroutine bc_finalize
460
463 subroutine bc_list_init(bclst, size)
464 type(bc_list_t), intent(inout), target :: bclst
465 integer, optional :: size
466 integer :: n, i
467
468 call bc_list_free(bclst)
469
470 if (present(size)) then
471 n = size
472 else
473 n = 1
474 end if
475
476 allocate(bclst%bc(n))
477
478 do i = 1, n
479 bclst%bc(i)%bcp => null()
480 end do
481
482 bclst%n = 0
483 bclst%size = n
484
485 end subroutine bc_list_init
486
490 subroutine bc_list_free(bclst)
491 type(bc_list_t), intent(inout) :: bclst
492
493 if (allocated(bclst%bc)) then
494 deallocate(bclst%bc)
495 end if
496
497 bclst%n = 0
498 bclst%size = 0
499
500 end subroutine bc_list_free
501
504 subroutine bc_list_add(bclst, bc)
505 type(bc_list_t), intent(inout) :: bclst
506 class(bc_t), intent(inout), target :: bc
507 type(bcp_t), allocatable :: tmp(:)
508
510 if(bc%marked_facet%size() .eq. 0) return
511
512 if (bclst%n .ge. bclst%size) then
513 bclst%size = bclst%size * 2
514 allocate(tmp(bclst%size))
515 tmp(1:bclst%n) = bclst%bc
516 call move_alloc(tmp, bclst%bc)
517 end if
518
519 bclst%n = bclst%n + 1
520 bclst%bc(bclst%n)%bcp => bc
521
522 end subroutine bc_list_add
523
529 subroutine bc_list_apply_scalar(bclst, x, n, t, tstep)
530 type(bc_list_t), intent(inout) :: bclst
531 integer, intent(in) :: n
532 real(kind=rp), intent(inout), dimension(n) :: x
533 real(kind=rp), intent(in), optional :: t
534 integer, intent(in), optional :: tstep
535 type(c_ptr) :: x_d
536 integer :: i
537
538 if (neko_bcknd_device .eq. 1) then
539 x_d = device_get_ptr(x)
540 if (present(t) .and. present(tstep)) then
541 do i = 1, bclst%n
542 call bclst%bc(i)%bcp%apply_scalar_dev(x_d, t=t, tstep=tstep)
543 end do
544 else if (present(t)) then
545 do i = 1, bclst%n
546 call bclst%bc(i)%bcp%apply_scalar_dev(x_d, t=t)
547 end do
548 else if (present(tstep)) then
549 do i = 1, bclst%n
550 call bclst%bc(i)%bcp%apply_scalar_dev(x_d, tstep=tstep)
551 end do
552 else
553 do i = 1, bclst%n
554 call bclst%bc(i)%bcp%apply_scalar_dev(x_d)
555 end do
556 end if
557 else
558 if (present(t) .and. present(tstep)) then
559 do i = 1, bclst%n
560 call bclst%bc(i)%bcp%apply_scalar(x, n, t, tstep)
561 end do
562 else if (present(t)) then
563 do i = 1, bclst%n
564 call bclst%bc(i)%bcp%apply_scalar(x, n, t=t)
565 end do
566 else if (present(tstep)) then
567 do i = 1, bclst%n
568 call bclst%bc(i)%bcp%apply_scalar(x, n, tstep=tstep)
569 end do
570 else
571 do i = 1, bclst%n
572 call bclst%bc(i)%bcp%apply_scalar(x, n)
573 end do
574 end if
575 end if
576
577 end subroutine bc_list_apply_scalar
578
586 subroutine bc_list_apply_vector(bclst, x, y, z, n, t, tstep)
587 type(bc_list_t), intent(inout) :: bclst
588 integer, intent(in) :: n
589 real(kind=rp), intent(inout), dimension(n) :: x
590 real(kind=rp), intent(inout), dimension(n) :: y
591 real(kind=rp), intent(inout), dimension(n) :: z
592 real(kind=rp), intent(in), optional :: t
593 integer, intent(in), optional :: tstep
594 type(c_ptr) :: x_d
595 type(c_ptr) :: y_d
596 type(c_ptr) :: z_d
597 integer :: i
598
599 if (neko_bcknd_device .eq. 1) then
600 x_d = device_get_ptr(x)
601 y_d = device_get_ptr(y)
602 z_d = device_get_ptr(z)
603 if (present(t) .and. present(tstep)) then
604 do i = 1, bclst%n
605 call bclst%bc(i)%bcp%apply_vector_dev(x_d, y_d, z_d, t, tstep)
606 end do
607 else if (present(t)) then
608 do i = 1, bclst%n
609 call bclst%bc(i)%bcp%apply_vector_dev(x_d, y_d, z_d, t=t)
610 end do
611 else if (present(tstep)) then
612 do i = 1, bclst%n
613 call bclst%bc(i)%bcp%apply_vector_dev(x_d, y_d, z_d, tstep=tstep)
614 end do
615 else
616 do i = 1, bclst%n
617 call bclst%bc(i)%bcp%apply_vector_dev(x_d, y_d, z_d)
618 end do
619 end if
620 else
621 if (present(t) .and. present(tstep)) then
622 do i = 1, bclst%n
623 call bclst%bc(i)%bcp%apply_vector(x, y, z, n, t, tstep)
624 end do
625 else if (present(t)) then
626 do i = 1, bclst%n
627 call bclst%bc(i)%bcp%apply_vector(x, y, z, n, t=t)
628 end do
629 else if (present(tstep)) then
630 do i = 1, bclst%n
631 call bclst%bc(i)%bcp%apply_vector(x, y, z, n, tstep=tstep)
632 end do
633 else
634 do i = 1, bclst%n
635 call bclst%bc(i)%bcp%apply_vector(x, y, z, n)
636 end do
637 end if
638 end if
639
640 end subroutine bc_list_apply_vector
641
642
643end module bc
Apply the boundary condition to a scalar field on the device.
Definition bc.f90:161
Apply the boundary condition to a scalar field.
Definition bc.f90:118
Apply the boundary condition to a vector field on the device.
Definition bc.f90:177
Apply the boundary condition to a vector field.
Definition bc.f90:137
Destructor.
Definition bc.f90:152
Return the device pointer for an associated Fortran array.
Definition device.F90:81
Map a Fortran array to a device (allocate and associate)
Definition device.F90:57
Copy data between host and device (or device and device)
Definition device.F90:51
Defines a boundary condition.
Definition bc.f90:34
subroutine bc_mark_zone(this, bc_zone)
Mark all facets from a zone.
Definition bc.f90:279
subroutine, public bc_list_add(bclst, bc)
Add a condition to a list of boundary conditions.
Definition bc.f90:505
subroutine, public bc_list_apply_vector(bclst, x, y, z, n, t, tstep)
Apply a list of boundary conditions to a vector field.
Definition bc.f90:587
subroutine bc_finalize(this)
Finalize the construction of the bc by populting the msk and facet arrays.
Definition bc.f90:364
subroutine bc_free_base(this)
Destructor for the base type, bc_t.
Definition bc.f90:218
subroutine, public bc_list_init(bclst, size)
Constructor for a list of boundary conditions.
Definition bc.f90:464
subroutine, public bc_list_free(bclst)
Destructor for a list of boundary conditions.
Definition bc.f90:491
subroutine bc_init_base(this, coef)
Constructor.
Definition bc.f90:202
subroutine bc_mark_zones_from_list(this, bc_zones, bc_key, bc_labels)
Mark all facets from a list of zones, also marks type of bc in the mesh. The facet_type in mesh is be...
Definition bc.f90:294
subroutine bc_mark_facet(this, facet, el)
Mark facet on element el as part of the boundary condition.
Definition bc.f90:251
subroutine bc_mark_facets(this, facet_list)
Mark all facets from a (facet, el) tuple list.
Definition bc.f90:264
subroutine, public bc_list_apply_scalar(bclst, x, n, t, tstep)
Apply a list of boundary conditions to a scalar field.
Definition bc.f90:530
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:185
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 mesh.
Definition mesh.f90:34
integer, parameter, public neko_msh_max_zlbls
Max num. zone labels.
Definition mesh.f90:56
integer, parameter, public neko_msh_max_zlbl_len
Max length of a zone label.
Definition mesh.f90:58
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
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:136
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:174
A list of boundary conditions.
Definition bc.f90:104
Base type for a boundary condition.
Definition bc.f90:51
Pointer to boundary condtiion.
Definition bc.f90:99
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
The function space for the SEM solution fields.
Definition space.f90:62
Integer 2-tuple based stack.
Definition stack.f90:84
Integer based 2-tuple.
Definition tuple.f90:51