Neko  0.9.99
A portable framework for high-order spectral element flow simulations
usr_scalar.f90
Go to the documentation of this file.
1 ! Copyright (c) 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 !
34 module usr_scalar
35  use num_types, only : rp
36  use coefs, only : coef_t
37  use bc, only : bc_t
38  use device
41  implicit none
42  private
43 
45  type, public, extends(bc_t) :: usr_scalar_t
46  procedure(usr_scalar_bc_eval), nopass, pointer :: eval => null()
47  type(c_ptr), private :: usr_x_d = c_null_ptr
48  contains
49  procedure, pass(this) :: apply_scalar => usr_scalar_apply_scalar
50  procedure, pass(this) :: apply_vector => usr_scalar_apply_vector
51  procedure, pass(this) :: validate => usr_scalar_validate
52  procedure, pass(this) :: set_eval => usr_scalar_set_eval
53  procedure, pass(this) :: apply_vector_dev => usr_scalar_apply_vector_dev
54  procedure, pass(this) :: apply_scalar_dev => usr_scalar_apply_scalar_dev
56  procedure, pass(this) :: free => usr_scalar_free
57  end type usr_scalar_t
58 
59  abstract interface
60 
77  subroutine usr_scalar_bc_eval(s, x, y, z, nx, ny, nz, &
78  ix, iy, iz, ie, t, tstep)
79  import rp
80  real(kind=rp), intent(inout) :: s
81  real(kind=rp), intent(in) :: x
82  real(kind=rp), intent(in) :: y
83  real(kind=rp), intent(in) :: z
84  real(kind=rp), intent(in) :: nx
85  real(kind=rp), intent(in) :: ny
86  real(kind=rp), intent(in) :: nz
87  integer, intent(in) :: ix
88  integer, intent(in) :: iy
89  integer, intent(in) :: iz
90  integer, intent(in) :: ie
91  real(kind=rp), intent(in) :: t
92  integer, intent(in) :: tstep
93  end subroutine usr_scalar_bc_eval
94  end interface
95 
96 
97  public :: usr_scalar_bc_eval
98 
99 contains
100 
101  subroutine usr_scalar_free(this)
102  class(usr_scalar_t), target, intent(inout) :: this
103 
104  call this%free_base
105 
106  if (c_associated(this%usr_x_d)) then
107  call device_free(this%usr_x_d)
108  end if
109 
110  end subroutine usr_scalar_free
111 
117  subroutine usr_scalar_apply_scalar(this, x, n, t, tstep)
118  class(usr_scalar_t), intent(inout) :: this
119  integer, intent(in) :: n
120  real(kind=rp), intent(inout), dimension(n) :: x
121  real(kind=rp), intent(in), optional :: t
122  integer, intent(in), optional :: tstep
123  integer :: i, m, k, idx(4), facet, tstep_
124  real(kind=rp) :: t_
125 
126  if (present(t)) then
127  t_ = t
128  else
129  t_ = 0.0_rp
130  end if
131 
132  if (present(tstep)) then
133  tstep_ = tstep
134  else
135  tstep_ = 1
136  end if
137 
138  associate(xc => this%coef%dof%x, yc => this%coef%dof%y, &
139  zc => this%coef%dof%z, nx => this%coef%nx, ny => this%coef%ny, &
140  nz => this%coef%nz, lx => this%coef%Xh%lx)
141  m = this%msk(0)
142  do i = 1, m
143  k = this%msk(i)
144  facet = this%facet(i)
145  idx = nonlinear_index(k, lx, lx, lx)
146  select case(facet)
147  case(1,2)
148  call this%eval(x(k), &
149  xc(idx(1), idx(2), idx(3), idx(4)), &
150  yc(idx(1), idx(2), idx(3), idx(4)), &
151  zc(idx(1), idx(2), idx(3), idx(4)), &
152  nx(idx(2), idx(3), facet, idx(4)), &
153  ny(idx(2), idx(3), facet, idx(4)), &
154  nz(idx(2), idx(3), facet, idx(4)), &
155  idx(1), idx(2), idx(3), idx(4), &
156  t_, tstep_)
157  case(3,4)
158  call this%eval(x(k), &
159  xc(idx(1), idx(2), idx(3), idx(4)), &
160  yc(idx(1), idx(2), idx(3), idx(4)), &
161  zc(idx(1), idx(2), idx(3), idx(4)), &
162  nx(idx(1), idx(3), facet, idx(4)), &
163  ny(idx(1), idx(3), facet, idx(4)), &
164  nz(idx(1), idx(3), facet, idx(4)), &
165  idx(1), idx(2), idx(3), idx(4), &
166  t_, tstep_)
167  case(5,6)
168  call this%eval(x(k), &
169  xc(idx(1), idx(2), idx(3), idx(4)), &
170  yc(idx(1), idx(2), idx(3), idx(4)), &
171  zc(idx(1), idx(2), idx(3), idx(4)), &
172  nx(idx(1), idx(2), facet, idx(4)), &
173  ny(idx(1), idx(2), facet, idx(4)), &
174  nz(idx(1), idx(2), facet, idx(4)), &
175  idx(1), idx(2), idx(3), idx(4), &
176  t_, tstep_)
177  end select
178  end do
179  end associate
180  end subroutine usr_scalar_apply_scalar
181 
187  subroutine usr_scalar_apply_scalar_dev(this, x_d, t, tstep)
188  class(usr_scalar_t), intent(inout), target :: this
189  type(c_ptr) :: x_d
190  real(kind=rp), intent(in), optional :: t
191  integer, intent(in), optional :: tstep
192  integer :: i, m, k, idx(4), facet, tstep_
193  real(kind=rp) :: t_
194  integer(c_size_t) :: s
195  real(kind=rp), allocatable :: x(:)
196  m = this%msk(0)
197 
198  if (present(t)) then
199  t_ = t
200  else
201  t_ = 0.0_rp
202  end if
203 
204  if (present(tstep)) then
205  tstep_ = tstep
206  else
207  tstep_ = 1
208  end if
209 
210  associate(xc => this%coef%dof%x, yc => this%coef%dof%y, &
211  zc => this%coef%dof%z, nx => this%coef%nx, ny => this%coef%ny, &
212  nz => this%coef%nz, lx => this%coef%Xh%lx, usr_x_d => this%usr_x_d)
213 
214 
215  ! Pretabulate values during first call to apply
216  if (.not. c_associated(usr_x_d)) then
217  allocate(x(m)) ! Temp arrays
218  s = m*rp
219 
220  call device_alloc(this%usr_x_d, s)
221 
222  do i = 1, m
223  k = this%msk(i)
224  facet = this%facet(i)
225  idx = nonlinear_index(k, lx, lx, lx)
226  select case(facet)
227  case(1,2)
228  call this%eval(x(i), &
229  xc(idx(1), idx(2), idx(3), idx(4)), &
230  yc(idx(1), idx(2), idx(3), idx(4)), &
231  zc(idx(1), idx(2), idx(3), idx(4)), &
232  nx(idx(2), idx(3), facet, idx(4)), &
233  ny(idx(2), idx(3), facet, idx(4)), &
234  nz(idx(2), idx(3), facet, idx(4)), &
235  idx(1), idx(2), idx(3), idx(4), &
236  t_, tstep_)
237  case(3,4)
238  call this%eval(x(i), &
239  xc(idx(1), idx(2), idx(3), idx(4)), &
240  yc(idx(1), idx(2), idx(3), idx(4)), &
241  zc(idx(1), idx(2), idx(3), idx(4)), &
242  nx(idx(1), idx(3), facet, idx(4)), &
243  ny(idx(1), idx(3), facet, idx(4)), &
244  nz(idx(1), idx(3), facet, idx(4)), &
245  idx(1), idx(2), idx(3), idx(4), &
246  t_, tstep_)
247  case(5,6)
248  call this%eval(x(i), &
249  xc(idx(1), idx(2), idx(3), idx(4)), &
250  yc(idx(1), idx(2), idx(3), idx(4)), &
251  zc(idx(1), idx(2), idx(3), idx(4)), &
252  nx(idx(1), idx(2), facet, idx(4)), &
253  ny(idx(1), idx(2), facet, idx(4)), &
254  nz(idx(1), idx(2), facet, idx(4)), &
255  idx(1), idx(2), idx(3), idx(4), &
256  t_, tstep_)
257  end select
258  end do
259 
260  call device_memcpy(x, this%usr_x_d, m, host_to_device, sync=.true.)
261 
262  deallocate(x)
263  end if
264 
265  call device_inhom_dirichlet_apply_scalar(this%msk_d, x_d, &
266  this%usr_x_d, m)
267  end associate
268 
269 
270  end subroutine usr_scalar_apply_scalar_dev
271 
273  subroutine usr_scalar_apply_vector(this, x, y, z, n, t, tstep)
274  class(usr_scalar_t), intent(inout) :: this
275  integer, intent(in) :: n
276  real(kind=rp), intent(inout), dimension(n) :: x
277  real(kind=rp), intent(inout), dimension(n) :: y
278  real(kind=rp), intent(inout), dimension(n) :: z
279  real(kind=rp), intent(in), optional :: t
280  integer, intent(in), optional :: tstep
281 
282  end subroutine usr_scalar_apply_vector
283 
285  subroutine usr_scalar_apply_vector_dev(this, x_d, y_d, z_d, t, tstep)
286  class(usr_scalar_t), intent(inout), target :: this
287  type(c_ptr) :: x_d
288  type(c_ptr) :: y_d
289  type(c_ptr) :: z_d
290  real(kind=rp), intent(in), optional :: t
291  integer, intent(in), optional :: tstep
292 
293  end subroutine usr_scalar_apply_vector_dev
294 
297  subroutine usr_scalar_set_eval(this, user_scalar_bc)
298  class(usr_scalar_t), intent(inout) :: this
299  procedure(usr_scalar_bc_eval) :: user_scalar_bc
300  this%eval => user_scalar_bc
301  end subroutine usr_scalar_set_eval
302 
304  subroutine usr_scalar_validate(this)
305  class(usr_scalar_t), intent(inout) :: this
306  logical :: valid
307 
308  valid = .true. ! Assert it's going to be ok...
309  if (.not. associated(this%coef)) then
310  call neko_warning('Missing coefficients')
311  valid = .false.
312  end if
313 
314  if (.not. associated(this%eval)) then
315  call neko_warning('Missing eval function')
316  valid = .false.
317  end if
318 
319  if (.not. valid) then
320  call neko_error('Invalid user defined scalar condition')
321  end if
322 
323  end subroutine usr_scalar_validate
324 
325 end module usr_scalar
__device__ void nonlinear_index(const int idx, const int lx, int *index)
Abstract interface defining a user defined scalar boundary condition (pointwise) Just imitating inflo...
Definition: usr_scalar.f90:77
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
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Defines dirichlet conditions for scalars.
Definition: usr_scalar.f90:34
subroutine usr_scalar_apply_scalar_dev(this, x_d, t, tstep)
Scalar apply (device version) Just imitating inflow for now, but we should look this over Applies bou...
Definition: usr_scalar.f90:188
subroutine usr_scalar_apply_vector_dev(this, x_d, y_d, z_d, t, tstep)
No-op vector apply (device version)
Definition: usr_scalar.f90:286
subroutine usr_scalar_apply_vector(this, x, y, z, n, t, tstep)
No-op vector apply.
Definition: usr_scalar.f90:274
subroutine usr_scalar_set_eval(this, user_scalar_bc)
Assign user provided eval function.
Definition: usr_scalar.f90:298
subroutine usr_scalar_apply_scalar(this, x, n, t, tstep)
Scalar apply Just imitating inflow for now, but we should look this over Applies boundary conditions ...
Definition: usr_scalar.f90:118
subroutine usr_scalar_validate(this)
Validate user scalar condition.
Definition: usr_scalar.f90:305
subroutine usr_scalar_free(this)
Definition: usr_scalar.f90:102
Utilities.
Definition: utils.f90:35
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Definition: utils.f90:245
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
User defined dirichlet condition for scalars.
Definition: usr_scalar.f90:45