Neko 0.9.99
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
usr_inflow.f90
Go to the documentation of this file.
1! Copyright (c) 2020-2023, 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!
35 use num_types, only : rp
36 use coefs, only : coef_t
37 use inflow, only : inflow_t
38 use device
41 use bc, only : bc_t
42 implicit none
43 private
44
46 type, public, extends(bc_t) :: usr_inflow_t
47 procedure(usr_inflow_eval), nopass, pointer :: eval => null()
48 type(c_ptr), private :: usr_x_d = c_null_ptr
49 type(c_ptr), private :: usr_y_d = c_null_ptr
50 type(c_ptr), private :: usr_z_d = c_null_ptr
51 contains
52 procedure, pass(this) :: apply_scalar => usr_inflow_apply_scalar
53 procedure, pass(this) :: apply_vector => usr_inflow_apply_vector
54 procedure, pass(this) :: validate => usr_inflow_validate
55 procedure, pass(this) :: set_eval => usr_inflow_set_eval
56 procedure, pass(this) :: apply_vector_dev => usr_inflow_apply_vector_dev
57 procedure, pass(this) :: apply_scalar_dev => usr_inflow_apply_scalar_dev
59 procedure, pass(this) :: free => usr_inflow_free
60 end type usr_inflow_t
61
62 abstract interface
63
80 subroutine usr_inflow_eval(u, v, w, x, y, z, nx, ny, nz, &
81 ix, iy, iz, ie, t, tstep)
82 import rp
83 real(kind=rp), intent(inout) :: u
84 real(kind=rp), intent(inout) :: v
85 real(kind=rp), intent(inout) :: w
86 real(kind=rp), intent(in) :: x
87 real(kind=rp), intent(in) :: y
88 real(kind=rp), intent(in) :: z
89 real(kind=rp), intent(in) :: nx
90 real(kind=rp), intent(in) :: ny
91 real(kind=rp), intent(in) :: nz
92 integer, intent(in) :: ix
93 integer, intent(in) :: iy
94 integer, intent(in) :: iz
95 integer, intent(in) :: ie
96 real(kind=rp), intent(in) :: t
97 integer, intent(in) :: tstep
98 end subroutine usr_inflow_eval
99 end interface
100
101 public :: usr_inflow_eval
102
103contains
104
105 subroutine usr_inflow_free(this)
106 class(usr_inflow_t), target, intent(inout) :: this
107
108 call this%free_base()
109
110 if (c_associated(this%usr_x_d)) then
111 call device_free(this%usr_x_d)
112 end if
113
114 if (c_associated(this%usr_y_d)) then
115 call device_free(this%usr_y_d)
116 end if
117
118 if (c_associated(this%usr_z_d)) then
119 call device_free(this%usr_z_d)
120 end if
121
122 end subroutine usr_inflow_free
123
125 subroutine usr_inflow_apply_scalar(this, x, n, t, tstep)
126 class(usr_inflow_t), intent(inout) :: this
127 integer, intent(in) :: n
128 real(kind=rp), intent(inout), dimension(n) :: x
129 real(kind=rp), intent(in), optional :: t
130 integer, intent(in), optional :: tstep
131 end subroutine usr_inflow_apply_scalar
132
134 subroutine usr_inflow_apply_scalar_dev(this, x_d, t, tstep)
135 class(usr_inflow_t), intent(inout), target :: this
136 type(c_ptr) :: x_d
137 real(kind=rp), intent(in), optional :: t
138 integer, intent(in), optional :: tstep
139 end subroutine usr_inflow_apply_scalar_dev
140
142 subroutine usr_inflow_apply_vector(this, x, y, z, n, t, tstep)
143 class(usr_inflow_t), intent(inout) :: this
144 integer, intent(in) :: n
145 real(kind=rp), intent(inout), dimension(n) :: x
146 real(kind=rp), intent(inout), dimension(n) :: y
147 real(kind=rp), intent(inout), dimension(n) :: z
148 real(kind=rp), intent(in), optional :: t
149 integer, intent(in), optional :: tstep
150 integer :: i, m, k, idx(4), facet, tstep_
151 real(kind=rp) :: t_
152
153 if (present(t)) then
154 t_ = t
155 else
156 t_ = 0.0_rp
157 end if
158
159 if (present(tstep)) then
160 tstep_ = tstep
161 else
162 tstep_ = 1
163 end if
164
165 associate(&
166 xc => this%coef%dof%x, yc => this%coef%dof%y, zc => this%coef%dof%z, &
167 nx => this%coef%nx, ny => this%coef%ny, nz => this%coef%nz, &
168 lx => this%coef%Xh%lx)
169 m = this%msk(0)
170 do i = 1, m
171 k = this%msk(i)
172 facet = this%facet(i)
173 idx = nonlinear_index(k, lx, lx, lx)
174 select case (facet)
175 case (1,2)
176 call this%eval(x(k), y(k), z(k), &
177 xc(idx(1), idx(2), idx(3), idx(4)), &
178 yc(idx(1), idx(2), idx(3), idx(4)), &
179 zc(idx(1), idx(2), idx(3), idx(4)), &
180 nx(idx(2), idx(3), facet, idx(4)), &
181 ny(idx(2), idx(3), facet, idx(4)), &
182 nz(idx(2), idx(3), facet, idx(4)), &
183 idx(1), idx(2), idx(3), idx(4), &
184 t_, tstep_)
185 case (3,4)
186 call this%eval(x(k), y(k), z(k), &
187 xc(idx(1), idx(2), idx(3), idx(4)), &
188 yc(idx(1), idx(2), idx(3), idx(4)), &
189 zc(idx(1), idx(2), idx(3), idx(4)), &
190 nx(idx(1), idx(3), facet, idx(4)), &
191 ny(idx(1), idx(3), facet, idx(4)), &
192 nz(idx(1), idx(3), facet, idx(4)), &
193 idx(1), idx(2), idx(3), idx(4), &
194 t_, tstep_)
195 case (5,6)
196 call this%eval(x(k), y(k), z(k), &
197 xc(idx(1), idx(2), idx(3), idx(4)), &
198 yc(idx(1), idx(2), idx(3), idx(4)), &
199 zc(idx(1), idx(2), idx(3), idx(4)), &
200 nx(idx(1), idx(2), facet, idx(4)), &
201 ny(idx(1), idx(2), facet, idx(4)), &
202 nz(idx(1), idx(2), facet, idx(4)), &
203 idx(1), idx(2), idx(3), idx(4), &
204 t_, tstep_)
205 end select
206 end do
207 end associate
208
209 end subroutine usr_inflow_apply_vector
210
211 subroutine usr_inflow_apply_vector_dev(this, x_d, y_d, z_d, t, tstep)
212 class(usr_inflow_t), intent(inout), target :: this
213 type(c_ptr) :: x_d
214 type(c_ptr) :: y_d
215 type(c_ptr) :: z_d
216 real(kind=rp), intent(in), optional :: t
217 integer, intent(in), optional :: tstep
218 integer :: i, m, k, idx(4), facet, tstep_
219 integer(c_size_t) :: s
220 real(kind=rp) :: t_
221 real(kind=rp), allocatable :: x(:)
222 real(kind=rp), allocatable :: y(:)
223 real(kind=rp), allocatable :: z(:)
224
225 if (present(t)) then
226 t_ = t
227 else
228 t_ = 0.0_rp
229 end if
230
231 if (present(tstep)) then
232 tstep_ = tstep
233 else
234 tstep_ = 1
235 end if
236
237 associate(&
238 xc => this%coef%dof%x, yc => this%coef%dof%y, zc => this%coef%dof%z, &
239 nx => this%coef%nx, ny => this%coef%ny, nz => this%coef%nz, &
240 lx => this%coef%Xh%lx, usr_x_d => this%usr_x_d, &
241 usr_y_d => this%usr_y_d, usr_z_d => this%usr_z_d)
242
243 m = this%msk(0)
244
245
246 ! Pretabulate values during first call to apply
247 if (.not. c_associated(usr_x_d)) then
248 allocate(x(m), y(m), z(m)) ! Temp arrays
249
250 s = m*rp
251
252 call device_alloc(usr_x_d, s)
253 call device_alloc(usr_y_d, s)
254 call device_alloc(usr_z_d, s)
255
256 associate(xc => this%coef%dof%x, yc => this%coef%dof%y, &
257 zc => this%coef%dof%z, &
258 nx => this%coef%nx, ny => this%coef%ny, nz => this%coef%nz, &
259 lx => this%coef%Xh%lx)
260 do i = 1, m
261 k = this%msk(i)
262 facet = this%facet(i)
263 idx = nonlinear_index(k, lx, lx, lx)
264 select case (facet)
265 case (1,2)
266 call this%eval(x(i), y(i), z(i), &
267 xc(idx(1), idx(2), idx(3), idx(4)), &
268 yc(idx(1), idx(2), idx(3), idx(4)), &
269 zc(idx(1), idx(2), idx(3), idx(4)), &
270 nx(idx(2), idx(3), facet, idx(4)), &
271 ny(idx(2), idx(3), facet, idx(4)), &
272 nz(idx(2), idx(3), facet, idx(4)), &
273 idx(1), idx(2), idx(3), idx(4), &
274 t_, tstep_)
275 case (3,4)
276 call this%eval(x(i), y(i), z(i), &
277 xc(idx(1), idx(2), idx(3), idx(4)), &
278 yc(idx(1), idx(2), idx(3), idx(4)), &
279 zc(idx(1), idx(2), idx(3), idx(4)), &
280 nx(idx(1), idx(3), facet, idx(4)), &
281 ny(idx(1), idx(3), facet, idx(4)), &
282 nz(idx(1), idx(3), facet, idx(4)), &
283 idx(1), idx(2), idx(3), idx(4), &
284 t_, tstep_)
285 case (5,6)
286 call this%eval(x(i), y(i), z(i), &
287 xc(idx(1), idx(2), idx(3), idx(4)), &
288 yc(idx(1), idx(2), idx(3), idx(4)), &
289 zc(idx(1), idx(2), idx(3), idx(4)), &
290 nx(idx(1), idx(2), facet, idx(4)), &
291 ny(idx(1), idx(2), facet, idx(4)), &
292 nz(idx(1), idx(2), facet, idx(4)), &
293 idx(1), idx(2), idx(3), idx(4), &
294 t_, tstep_)
295 end select
296 end do
297 end associate
298
299 call device_memcpy(x, usr_x_d, m, host_to_device, sync = .false.)
300 call device_memcpy(y, usr_y_d, m, host_to_device, sync = .false.)
301 call device_memcpy(z, usr_z_d, m, host_to_device, sync = .true.)
302
303 deallocate(x, y, z)
304 end if
305
306 call device_inhom_dirichlet_apply_vector(this%msk_d, x_d, y_d, z_d, &
307 usr_x_d, usr_y_d, usr_z_d, m)
308
309 end associate
310
311 end subroutine usr_inflow_apply_vector_dev
312
315 subroutine usr_inflow_set_eval(this, usr_eval)
316 class(usr_inflow_t), intent(inout) :: this
317 procedure(usr_inflow_eval) :: usr_eval
318 this%eval => usr_eval
319 end subroutine usr_inflow_set_eval
320
322 subroutine usr_inflow_validate(this)
323 class(usr_inflow_t), intent(inout) :: this
324 logical :: valid
325
326 valid = .true. ! Assert it's going to be ok...
327 if (.not. associated(this%coef)) then
328 call neko_warning('Missing coefficients')
329 valid = .false.
330 end if
331
332 if (.not. associated(this%eval)) then
333 call neko_warning('Missing eval function')
334 valid = .false.
335 end if
336
337 if (.not. valid) then
338 call neko_error('Invalid user defined inflow condition')
339 end if
340
341 end subroutine usr_inflow_validate
342
343end module usr_inflow
__device__ void nonlinear_index(const int idx, const int lx, int *index)
Abstract interface defining a user defined inflow condition (pointwise)
Defines a boundary condition.
Definition bc.f90:34
Coefficients.
Definition coef.f90:34
Device abstraction, common interface for various accelerators.
Definition device.F90:34
subroutine, public device_free(x_d)
Deallocate memory on the device.
Definition device.F90:185
Defines inflow dirichlet conditions.
Definition inflow.f90:34
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Defines inflow dirichlet conditions.
subroutine usr_inflow_apply_scalar(this, x, n, t, tstep)
No-op scalar apply.
subroutine usr_inflow_apply_vector_dev(this, x_d, y_d, z_d, t, tstep)
subroutine usr_inflow_apply_vector(this, x, y, z, n, t, tstep)
Apply user defined inflow conditions (vector valued)
subroutine usr_inflow_apply_scalar_dev(this, x_d, t, tstep)
No-op scalar apply (device version)
subroutine usr_inflow_validate(this)
Validate user inflow condition.
subroutine usr_inflow_free(this)
subroutine usr_inflow_set_eval(this, usr_eval)
Assign user provided eval function.
Utilities.
Definition utils.f90:35
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Definition utils.f90:266
Base type for a boundary condition.
Definition bc.f90:51
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
Dirichlet condition for inlet (vector valued)
Definition inflow.f90:43
User defined dirichlet condition for inlet (vector valued)