Neko 0.9.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
gradient_jump_penalty.f90
Go to the documentation of this file.
1! Copyright (c) 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!
33!
35module gradient_jump_penalty
36 use num_types, only : rp
37 use utils, only : neko_error
38 use json_utils, only : json_get_or_default
39 use json_module, only : json_file
40 use math
41 use point, only : point_t
42 use field, only : field_t
43 use dofmap , only : dofmap_t
44 use neko_config, only : neko_bcknd_device
45 use coefs, only : coef_t
46 use element, only : element_t
47 use hex
48 use quad
49 use operators, only : dudxyz
50 use gs_ops, only : gs_op_add
51 use space
52 use gather_scatter, only : gs_t
53 use device
54 use device_math
55 use device_gradient_jump_penalty
56
57 implicit none
58 private
59
62 type, public :: gradient_jump_penalty_t
64 real(kind=rp) :: tau
66 integer :: p
68 integer :: lx
70 real(kind=rp), allocatable, dimension(:, :, :, :) :: penalty
71 type(c_ptr) :: penalty_d = c_null_ptr
73 real(kind=rp), allocatable, dimension(:, :, :, :) :: penalty_facet
74 type(c_ptr) :: penalty_facet_d = c_null_ptr
76 type(coef_t), pointer :: coef => null()
78 type(dofmap_t), pointer :: dm => null()
80 real(kind=rp), allocatable, dimension(:, :, :, :) :: grad1, &
81 grad2, grad3
82 type(c_ptr) :: grad1_d = c_null_ptr
83 type(c_ptr) :: grad2_d = c_null_ptr
84 type(c_ptr) :: grad3_d = c_null_ptr
86 real(kind=rp), allocatable, dimension(:, :, :, :) :: g
87 type(c_ptr) :: g_d = c_null_ptr
89 real(kind=rp), allocatable, dimension(:, :, :, :) :: flux1, flux2, flux3
90 type(c_ptr) :: flux1_d = c_null_ptr
91 type(c_ptr) :: flux2_d = c_null_ptr
92 type(c_ptr) :: flux3_d = c_null_ptr
94 real(kind=rp), allocatable, dimension(:, :, :, :) :: volflux1, &
95 volflux2, volflux3
96 type(c_ptr) :: volflux1_d = c_null_ptr
97 type(c_ptr) :: volflux2_d = c_null_ptr
98 type(c_ptr) :: volflux3_d = c_null_ptr
100 real(kind=rp), allocatable, dimension(:, :, :, :) :: absvolflux
101 type(c_ptr) :: absvolflux_d = c_null_ptr
103 real(kind=rp), allocatable, dimension(:, :, :, :) :: n1, n2, n3
104 type(c_ptr) :: n1_d = c_null_ptr
105 type(c_ptr) :: n2_d = c_null_ptr
106 type(c_ptr) :: n3_d = c_null_ptr
108 real(kind=rp), allocatable, dimension(:, :, :, :) :: facet_factor
109 type(c_ptr) :: facet_factor_d = c_null_ptr
111 integer, allocatable :: n_facet(:)
112 integer :: n_facet_max
114 real(kind=rp), allocatable, dimension(:, :, :, :) :: h2
116 real(kind=rp), allocatable :: dphidxi(:, :)
117 type(c_ptr) :: dphidxi_d = c_null_ptr
119 type(space_t) :: xh_gjp
120 type(dofmap_t) :: dm_gjp
121 type(gs_t) :: gs_gjp
123 integer :: n
125 integer :: n_large
126
127 contains
129 procedure, pass(this) :: init => gradient_jump_penalty_init
131 procedure, pass(this) :: free => gradient_jump_penalty_free
133 procedure, pass(this) :: compute => gradient_jump_penalty_compute
135 procedure, pass(this) :: perform => gradient_jump_penalty_perform
136
137 end type gradient_jump_penalty_t
138
139contains
145 subroutine gradient_jump_penalty_init(this, params, dofmap, coef, a, b)
146 implicit none
147 class(gradient_jump_penalty_t), intent(inout) :: this
148 type(json_file), target, intent(inout) :: params
149 type(dofmap_t), target, intent(in) :: dofmap
150 type(coef_t), target, intent(in) :: coef
151 real(kind=rp), intent(in) :: a, b
152
153 integer :: i, j, k, l
154 real(kind=rp), allocatable :: zg(:) ! Quadrature points
155 real(kind=rp) :: normal(3)
156
157 call this%free()
158
159 this%dm => dofmap
160 this%p = dofmap%xh%lx - 1
161 this%lx = dofmap%xh%lx
162
163 if (this%p .gt. 1) then
164 this%tau = -a * (this%p + 1) ** (-b)
165 else
166 this%tau = -a
167 end if
168
169 this%coef => coef
170 this%n = this%lx ** 3 * this%coef%msh%nelv
171 this%n_large = (this%lx + 2) ** 3 * this%coef%msh%nelv
172
173 allocate(this%n_facet(this%coef%msh%nelv))
174 do i = 1, this%coef%msh%nelv
175 select type (ep => this%coef%msh%elements(i)%e)
176 type is (hex_t)
177 this%n_facet(i) = 6
178 type is (quad_t)
179 call neko_error("Only Hexahedral element is &
180 &supported now for gradient jump penalty")
181 end select
182 end do
183 this%n_facet_max = maxval(this%n_facet)
184
185 allocate(this%h2(this%lx + 2, this%lx + 2, &
186 this%lx + 2, this%coef%msh%nelv))
187
188 do i = 1, this%coef%msh%nelv
189 select type (ep => this%coef%msh%elements(i)%e)
190 type is (hex_t)
191 call eval_h2_hex(this%h2(:, :, :, i), this%lx, i, this%dm, this%coef)
192 type is (quad_t)
193 call neko_error("Gradient jump penalty error: mesh size &
194 &evaluation is not supported for quad_t")
195 end select
196 end do
197
198 allocate(zg(this%lx))
199 allocate(this%dphidxi(this%lx, this%lx))
200
201 zg = dofmap%xh%zg(:,1)
202 do i = 1, dofmap%xh%lx
203 do j = 1, dofmap%xh%lx
204 this%dphidxi(j,i) = this%coef%Xh%dx(j,i)
205 end do
206 end do
207
208 allocate(this%penalty(this%lx, this%lx, this%lx, this%coef%msh%nelv))
209 allocate(this%grad1(this%lx, this%lx, this%lx, this%coef%msh%nelv))
210 allocate(this%grad2(this%lx, this%lx, this%lx, this%coef%msh%nelv))
211 allocate(this%grad3(this%lx, this%lx, this%lx, this%coef%msh%nelv))
212
213 allocate(this%penalty_facet(this%lx + 2, this%lx + 2, &
214 this%lx + 2, this%coef%msh%nelv))
215 allocate(this%G(this%lx + 2, this%lx + 2, &
216 this%lx + 2, this%coef%msh%nelv))
217 allocate(this%flux1(this%lx + 2, this%lx + 2, &
218 this%lx + 2, this%coef%msh%nelv))
219 allocate(this%flux2(this%lx + 2, this%lx + 2, &
220 this%lx + 2, this%coef%msh%nelv))
221 allocate(this%flux3(this%lx + 2, this%lx + 2, &
222 this%lx + 2, this%coef%msh%nelv))
223 allocate(this%volflux1(this%lx + 2, this%lx + 2, &
224 this%lx + 2, this%coef%msh%nelv))
225 allocate(this%volflux2(this%lx + 2, this%lx + 2, &
226 this%lx + 2, this%coef%msh%nelv))
227 allocate(this%volflux3(this%lx + 2, this%lx + 2, &
228 this%lx + 2, this%coef%msh%nelv))
229 allocate(this%absvolflux(this%lx + 2, this%lx + 2, &
230 this%lx + 2, this%coef%msh%nelv))
231 allocate(this%n1(this%lx + 2, this%lx + 2, &
232 this%lx + 2, this%coef%msh%nelv))
233 allocate(this%n2(this%lx + 2, this%lx + 2, &
234 this%lx + 2, this%coef%msh%nelv))
235 allocate(this%n3(this%lx + 2, this%lx + 2, &
236 this%lx + 2, this%coef%msh%nelv))
237
238 ! Extract facets' normals
239 do i = 1, this%coef%msh%nelv
240 do j = 1, 6 ! for hexahedral elements
241 do k = 1, this%lx
242 do l = 1, this%lx
243 select case(j)
244 case(1)
245 normal = this%coef%get_normal(1, l, k, i, j)
246 this%n1(1, l + 1, k + 1, i) = normal(1)
247 this%n2(1, l + 1, k + 1, i) = normal(2)
248 this%n3(1, l + 1, k + 1, i) = normal(3)
249 case(2)
250 normal = this%coef%get_normal(1, l, k, i, j)
251 this%n1(this%lx + 2, l + 1, k + 1, i) = normal(1)
252 this%n2(this%lx + 2, l + 1, k + 1, i) = normal(2)
253 this%n3(this%lx + 2, l + 1, k + 1, i) = normal(3)
254 case(3)
255 normal = this%coef%get_normal(l, 1, k, i, j)
256 this%n1(l + 1, 1, k + 1, i) = normal(1)
257 this%n2(l + 1, 1, k + 1, i) = normal(2)
258 this%n3(l + 1, 1, k + 1, i) = normal(3)
259 case(4)
260 normal = this%coef%get_normal(l, 1, k, i, j)
261 this%n1(l + 1, this%lx + 2, k + 1, i) = normal(1)
262 this%n2(l + 1, this%lx + 2, k + 1, i) = normal(2)
263 this%n3(l + 1, this%lx + 2, k + 1, i) = normal(3)
264 case(5)
265 normal = this%coef%get_normal(l, k, 1, i, j)
266 this%n1(l + 1, k + 1, 1, i) = normal(1)
267 this%n2(l + 1, k + 1, 1, i) = normal(2)
268 this%n3(l + 1, k + 1, 1, i) = normal(3)
269 case(6)
270 normal = this%coef%get_normal(l, k, 1, i, j)
271 this%n1(l + 1, k + 1, this%lx + 2, i) = normal(1)
272 this%n2(l + 1, k + 1, this%lx + 2, i) = normal(2)
273 this%n3(l + 1, k + 1, this%lx + 2, i) = normal(3)
274 case default
275 call neko_error("The face index is not correct")
276 end select
277 end do
278 end do
279 end do
280 end do
281
282 ! Assemble facet factor
283 call facet_factor_init(this)
284
285 ! Initialize Gather-Scatter
286 call this%Xh_GJP%init(gll, this%lx+2, this%lx+2, this%lx+2)
287 call this%dm_GJP%init(this%coef%msh, this%Xh_GJP)
288 call this%gs_GJP%init(this%dm_GJP)
289
290 ! Initialize pointers for device
291 if (neko_bcknd_device .eq. 1) then
292 call device_map(this%dphidxi, this%dphidxi_d, &
293 this%lx * this%lx)
294 call device_map(this%penalty, this%penalty_d, this%n)
295 call device_map(this%grad1, this%grad1_d, this%n)
296 call device_map(this%grad2, this%grad2_d, this%n)
297 call device_map(this%grad3, this%grad3_d, this%n)
298
299 call device_map(this%penalty_facet, this%penalty_facet_d, this%n_large)
300 call device_map(this%G, this%G_d, this%n_large)
301 call device_map(this%flux1, this%flux1_d, this%n_large)
302 call device_map(this%flux2, this%flux2_d, this%n_large)
303 call device_map(this%flux3, this%flux3_d, this%n_large)
304
305 call device_map(this%volflux1, this%volflux1_d, this%n_large)
306 call device_map(this%volflux2, this%volflux2_d, this%n_large)
307 call device_map(this%volflux3, this%volflux3_d, this%n_large)
308 call device_map(this%absvolflux, this%absvolflux_d, this%n_large)
309
310 call device_map(this%n1, this%n1_d, this%n_large)
311 call device_map(this%n2, this%n2_d, this%n_large)
312 call device_map(this%n3, this%n3_d, this%n_large)
313 call device_map(this%facet_factor, this%facet_factor_d, this%n_large)
314
315 call device_memcpy(this%dphidxi, this%dphidxi_d, &
316 this%lx * this%lx, &
317 host_to_device, sync = .false.)
318 call device_memcpy(this%n1, this%n1_d, this%n_large, &
319 host_to_device, sync = .false.)
320 call device_memcpy(this%n2, this%n2_d, this%n_large, &
321 host_to_device, sync = .false.)
322 call device_memcpy(this%n3, this%n3_d, this%n_large, &
323 host_to_device, sync = .false.)
324 call device_memcpy(this%facet_factor, this%facet_factor_d, this%n_large,&
325 host_to_device, sync = .false.)
326
327 end if
328
329 end subroutine gradient_jump_penalty_init
330
334 subroutine eval_h2_hex(h2_el, n, i, dm, coef)
335 integer, intent(in) :: n, i
336 type(dofmap_t), pointer, intent(in) :: dm
337 type(coef_t), pointer, intent(in) :: coef
338 real(kind=rp), intent(inout) :: h2_el(n + 2, n + 2, n + 2)
339
340 integer :: j, k, l
341
342 h2_el = 0.0_rp
343
344 do j = 1, 6
345 do k = 1, n
346 do l = 1, n
347 select case(j)
348 case(1)
349 h2_el(1, l + 1, k + 1) = &
350 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
351 case(2)
352 h2_el(n + 2, l + 1, k + 1) = &
353 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
354 case(3)
355 h2_el(l + 1, 1, k + 1) = &
356 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
357 case(4)
358 h2_el(l + 1, n + 2, k + 1) = &
359 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
360 case(5)
361 h2_el(l + 1, k + 1, 1) = &
362 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
363 case(6)
364 h2_el(l + 1, k + 1, n + 2) = &
365 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
366 case default
367 call neko_error("The face index is not correct")
368 end select
369 end do
370 end do
371 end do
372
373 end subroutine eval_h2_hex
374
375 function dist2_quadrature_hex(l, k, j, i, n, dm, coef) result(dist2)
376 integer, intent(in) :: l, k, j, i, n
377 type(dofmap_t), pointer, intent(in) :: dm
378 type(coef_t), pointer, intent(in) :: coef
379 real(kind=rp) :: dist2, dist_1, dist_2
380
381 real(kind=rp) :: x1, y1, z1, x2, y2, z2
382 real(kind=rp) :: normal1(3), normal2(3)
383 real(kind=rp) :: n11, n12, n13, n21, n22, n23
384 real(kind=rp) :: v1, v2, v3
385
386 dist2 = 0.0_rp
387 select case(j)
388 case(1)
389 normal1 = coef%get_normal(1, l, k, i, 1)
390 n11 = normal1(1)
391 n12 = normal1(2)
392 n13 = normal1(3)
393 x1 = dm%x(1, l, k, i)
394 y1 = dm%y(1, l, k, i)
395 z1 = dm%z(1, l, k, i)
396 normal2 = coef%get_normal(1, l, k, i, 2)
397 n21 = normal2(1)
398 n22 = normal2(2)
399 n23 = normal2(3)
400 x2 = dm%x(n, l, k, i)
401 y2 = dm%y(n, l, k, i)
402 z2 = dm%z(n, l, k, i)
403 case(2)
404 ! now the facet pair share the same value for h
405 ! but just let it be here for furture possible changes
406 normal1 = coef%get_normal(1, l, k, i, 2)
407 n11 = normal1(1)
408 n12 = normal1(2)
409 n13 = normal1(3)
410 x1 = dm%x(n, l, k, i)
411 y1 = dm%y(n, l, k, i)
412 z1 = dm%z(n, l, k, i)
413 normal2 = coef%get_normal(1, l, k, i, 1)
414 n21 = normal2(1)
415 n22 = normal2(2)
416 n23 = normal2(3)
417 x2 = dm%x(1, l, k, i)
418 y2 = dm%y(1, l, k, i)
419 z2 = dm%z(1, l, k, i)
420 case(3)
421 normal1 = coef%get_normal(1, l, k, i, 3)
422 n11 = normal1(1)
423 n12 = normal1(2)
424 n13 = normal1(3)
425 x1 = dm%x(l, 1, k, i)
426 y1 = dm%y(l, 1, k, i)
427 z1 = dm%z(l, 1, k, i)
428 normal2 = coef%get_normal(1, l, k, i, 4)
429 n21 = normal2(1)
430 n22 = normal2(2)
431 n23 = normal2(3)
432 x2 = dm%x(l, n, k, i)
433 y2 = dm%y(l, n, k, i)
434 z2 = dm%z(l, n, k, i)
435 case(4)
436 normal1 = coef%get_normal(1, l, k, i, 4)
437 n11 = normal1(1)
438 n12 = normal1(2)
439 n13 = normal1(3)
440 x1 = dm%x(l, n, k, i)
441 y1 = dm%y(l, n, k, i)
442 z1 = dm%z(l, n, k, i)
443 normal2 = coef%get_normal(1, l, k, i, 3)
444 n21 = normal2(1)
445 n22 = normal2(2)
446 n23 = normal2(3)
447 x2 = dm%x(l, 1, k, i)
448 y2 = dm%y(l, 1, k, i)
449 z2 = dm%z(l, 1, k, i)
450 case(5)
451 normal1 = coef%get_normal(1, l, k, i, 5)
452 n11 = normal1(1)
453 n12 = normal1(2)
454 n13 = normal1(3)
455 x1 = dm%x(l, k, 1, i)
456 y1 = dm%y(l, k, 1, i)
457 z1 = dm%z(l, k, 1, i)
458 normal2 = coef%get_normal(1, l, k, i, 6)
459 n21 = normal2(1)
460 n22 = normal2(2)
461 n23 = normal2(3)
462 x2 = dm%x(l, k, n, i)
463 y2 = dm%y(l, k, n, i)
464 z2 = dm%z(l, k, n, i)
465 case(6)
466 normal1 = coef%get_normal(1, l, k, i, 6)
467 n11 = normal1(1)
468 n12 = normal1(2)
469 n13 = normal1(3)
470 x1 = dm%x(l, k, n, i)
471 y1 = dm%y(l, k, n, i)
472 z1 = dm%z(l, k, n, i)
473 normal2 = coef%get_normal(1, l, k, i, 5)
474 n21 = normal2(1)
475 n22 = normal2(2)
476 n23 = normal2(3)
477 x2 = dm%x(l, k, 1, i)
478 y2 = dm%y(l, k, 1, i)
479 z2 = dm%z(l, k, 1, i)
480 case default
481 call neko_error("The face index is not correct")
482 end select
483
484 ! get the vector from the quadrature point to the one on the other side
485 v1 = x2 - x1
486 v2 = y2 - y1
487 v3 = z2 - z1
488 ! Project onto tabsvolflhe facet-normal direction of the point
489 dist_1 = v1*n11 + v2*n12 + v3*n13
490 dist_2 = - (v1*n21 + v2*n22 + v3*n23)
491
492 dist2 = ((dist_1 + dist_2)/2.0_rp)*((dist_1 + dist_2)/2.0_rp)
493
494 end function
495
497 subroutine facet_factor_init(this)
498 class(gradient_jump_penalty_t), intent(inout) :: this
499 integer :: i, j, k, l
500 real(kind=rp) :: area_tmp
501
502 allocate(this%facet_factor(this%lx + 2, this%lx + 2, &
503 this%lx + 2, this%coef%msh%nelv))
504
505 associate(facet_factor => this%facet_factor, &
506 coef => this%coef, &
507 lx => this%lx, &
508 nelv => this%coef%msh%nelv, &
509 jacinv => this%coef%jacinv, h2 => this%h2, &
510 tau => this%tau, n1 => this%n1, &
511 n2 => this%n2, n3 => this%n3)
512
513 do i = 1, nelv
514 do j = 1, 6 ! for hexahedral elementsh2
515 do k = 1, lx
516 do l = 1, lx
517 select case(j)
518 case(1)
519 area_tmp = coef%get_area(1, l, k, i, j)
520 facet_factor(1, l + 1, k + 1, i) = area_tmp * tau * &
521 h2(1, l + 1, k + 1, i) * &
522 (n1(1, l + 1, k + 1, i) * coef%drdx(1, l, k, i) + &
523 n2(1, l + 1, k + 1, i) * coef%drdy(1, l, k, i) + &
524 n3(1, l + 1, k + 1, i) * coef%drdz(1, l, k, i) ) &
525 * jacinv(1, l, k, i)
526 case(2)
527 area_tmp = coef%get_area(1, l, k, i, j)
528 facet_factor(lx + 2, l + 1, k + 1, i) = area_tmp * tau * &
529 h2(lx + 2, l + 1, k + 1, i) * &
530 (n1(lx + 2, l + 1, k + 1, i) * &
531 coef%drdx(lx, l, k, i) + &
532 n2(lx + 2, l + 1, k + 1, i) * &
533 coef%drdy(lx, l, k, i) + &
534 n3(lx + 2, l + 1, k + 1, i) * &
535 coef%drdz(lx, l, k, i) ) &
536 * jacinv(lx, l, k, i)
537 case(3)
538 area_tmp = coef%get_area(l, 1, k, i, j)
539 facet_factor(l + 1, 1, k + 1, i) = area_tmp * tau * &
540 h2(l + 1, 1, k + 1, i) * &
541 (n1(l + 1, 1, k + 1, i) * coef%dsdx(l, 1, k, i) + &
542 n2(l + 1, 1, k + 1, i) * coef%dsdy(l, 1, k, i) + &
543 n3(l + 1, 1, k + 1, i) * coef%dsdz(l, 1, k, i) ) &
544 * jacinv(l, 1, k, i)
545 case(4)
546 area_tmp = coef%get_area(l, 1, k, i, j)
547 facet_factor(l + 1, lx + 2, k + 1, i) = area_tmp * tau * &
548 h2(l + 1, lx + 2, k + 1, i) * &
549 (n1(l + 1, lx + 2, k + 1, i) * &
550 coef%dsdx(l, lx, k, i) + &
551 n2(l + 1, lx + 2, k + 1, i) * &
552 coef%dsdy(l, lx, k, i) + &
553 n3(l + 1, lx + 2, k + 1, i) * &
554 coef%dsdz(l, lx, k, i) ) &
555 * jacinv(l, lx, k, i)
556 case(5)
557 area_tmp = coef%get_area(l, k, 1, i, j)
558 facet_factor(l + 1, k + 1, 1, i) = area_tmp * tau * &
559 h2(l + 1, k + 1, 1, i) * &
560 (n1(l + 1, k + 1, 1, i) * coef%dtdx(l, k, 1, i) + &
561 n2(l + 1, k + 1, 1, i) * coef%dtdy(l, k, 1, i) + &
562 n3(l + 1, k + 1, 1, i) * coef%dtdz(l, k, 1, i) ) &
563 * jacinv(l, k, 1, i)
564 case(6)
565 area_tmp = coef%get_area(l, k, 1, i, j)
566 facet_factor(l + 1, k + 1, lx + 2, i) = area_tmp * tau * &
567 h2(l + 1, k + 1, lx + 2, i) * &
568 (n1(l + 1, k + 1, lx + 2, i) * coef%dtdx(l, k, lx, i) + &
569 n2(l + 1, k + 1, lx + 2, i) * coef%dtdy(l, k, lx, i) + &
570 n3(l + 1, k + 1, lx + 2, i) * coef%dtdz(l, k, lx, i) ) &
571 * jacinv(l, k, lx, i)
572 case default
573 call neko_error("The face index is not correct")
574 end select
575 end do
576 end do
577 end do
578 end do
579
580 end associate
581 end subroutine facet_factor_init
582
584 subroutine gradient_jump_penalty_free(this)
585 implicit none
586 class(gradient_jump_penalty_t), intent(inout) :: this
587
588 if (c_associated(this%dphidxi_d)) then
589 call device_free(this%dphidxi_d)
590 end if
591 if (c_associated(this%penalty_d)) then
592 call device_free(this%penalty_d)
593 end if
594 if (c_associated(this%grad1_d)) then
595 call device_free(this%grad1_d)
596 end if
597 if (c_associated(this%grad2_d)) then
598 call device_free(this%grad2_d)
599 end if
600 if (c_associated(this%grad3_d)) then
601 call device_free(this%grad3_d)
602 end if
603 if (c_associated(this%penalty_facet_d)) then
604 call device_free(this%penalty_facet_d)
605 end if
606 if (c_associated(this%G_d)) then
607 call device_free(this%G_d)
608 end if
609 if (c_associated(this%flux1_d)) then
610 call device_free(this%flux1_d)
611 end if
612 if (c_associated(this%flux2_d)) then
613 call device_free(this%flux2_d)
614 end if
615 if (c_associated(this%flux3_d)) then
616 call device_free(this%flux3_d)
617 end if
618 if (c_associated(this%volflux1_d)) then
619 call device_free(this%volflux1_d)
620 end if
621 if (c_associated(this%volflux2_d)) then
622 call device_free(this%volflux2_d)
623 end if
624 if (c_associated(this%volflux3_d)) then
625 call device_free(this%volflux3_d)
626 end if
627 if (c_associated(this%absvolflux_d)) then
628 call device_free(this%absvolflux_d)
629 end if
630 if (c_associated(this%n1_d)) then
631 call device_free(this%n1_d)
632 end if
633 if (c_associated(this%n2_d)) then
634 call device_free(this%n2_d)
635 end if
636 if (c_associated(this%n3_d)) then
637 call device_free(this%n3_d)
638 end if
639 if (c_associated(this%facet_factor_d)) then
640 call device_free(this%facet_factor_d)
641 end if
642
643 if (allocated(this%penalty)) then
644 deallocate(this%penalty)
645 end if
646 if (allocated(this%grad1)) then
647 deallocate(this%grad1)
648 end if
649 if (allocated(this%grad2)) then
650 deallocate(this%grad2)
651 end if
652 if (allocated(this%grad3)) then
653 deallocate(this%grad3)
654 end if
655 if (allocated(this%h2)) then
656 deallocate(this%h2)
657 end if
658 if (allocated(this%n_facet)) then
659 deallocate(this%n_facet)
660 end if
661 if (allocated(this%dphidxi)) then
662 deallocate(this%dphidxi)
663 end if
664 if (allocated(this%penalty_facet)) then
665 deallocate(this%penalty_facet)
666 end if
667 if (allocated(this%G)) then
668 deallocate(this%G)
669 end if
670 if (allocated(this%flux1)) then
671 deallocate(this%flux1)
672 end if
673 if (allocated(this%flux2)) then
674 deallocate(this%flux2)
675 end if
676 if (allocated(this%flux3)) then
677 deallocate(this%flux3)
678 end if
679 if (allocated(this%volflux1)) then
680 deallocate(this%volflux1)
681 end if
682 if (allocated(this%volflux2)) then
683 deallocate(this%volflux2)
684 end if
685 if (allocated(this%volflux3)) then
686 deallocate(this%volflux3)
687 end if
688 if (allocated(this%absvolflux)) then
689 deallocate(this%absvolflux)
690 end if
691 if (allocated(this%n1)) then
692 deallocate(this%n1)
693 end if
694 if (allocated(this%n2)) then
695 deallocate(this%n2)
696 end if
697 if (allocated(this%n3)) then
698 deallocate(this%n3)
699 end if
700 if (allocated(this%facet_factor)) then
701 deallocate(this%facet_factor)
702 end if
703
704 nullify(this%coef)
705 nullify(this%dm)
706
707 call this%Xh_GJP%free()
708 call this%gs_GJP%free()
709
710 end subroutine gradient_jump_penalty_free
711
717 subroutine gradient_jump_penalty_compute(this, u, v, w, s)
718 class(gradient_jump_penalty_t), intent(inout) :: this
719 type(field_t), intent(in) :: u, v, w, s
720
721 class(element_t), pointer :: ep
722 integer :: i
723
724 call g_compute(this, s)
725 call absvolflux_compute(this, u, v, w)
726
727 if (neko_bcknd_device .eq. 1) then
728 call device_col3(this%penalty_facet_d, this%absvolflux_d, this%G_d, &
729 this%n_large)
730 call device_col2(this%penalty_facet_d, this%facet_factor_d, this%n_large)
731 call device_gradient_jump_penalty_finalize(this%penalty_d, &
732 this%penalty_facet_d, &
733 this%dphidxi_d, &
734 this%lx, this%coef%msh%nelv)
735 else
736 call col3(this%penalty_facet, this%absvolflux, this%G, this%n_large)
737 call col2(this%penalty_facet, this%facet_factor, this%n_large)
738 call gradient_jump_penalty_finalize(this%penalty, this%penalty_facet, &
739 this%dphidxi, &
740 this%lx, this%coef%msh%nelv)
741 end if
742
743 end subroutine gradient_jump_penalty_compute
744
747 subroutine gradient_jump_penalty_perform(this, f)
748 class(gradient_jump_penalty_t), intent(inout) :: this
749 type(field_t), intent(inout) :: f
750
751 integer :: i, j, k, l
752
753 if (neko_bcknd_device .eq. 1) then
754 call device_add2(f%x_d, this%penalty_d, this%coef%dof%size())
755 else
756 call add2(f%x, this%penalty, this%coef%dof%size())
757 end if
758
759 end subroutine gradient_jump_penalty_perform
760
768 subroutine gradient_jump_penalty_finalize(penalty, wa, dphidxi, lx, nelv)
769 integer, intent(in) :: lx, nelv
770 real(kind=rp), intent(inout) :: penalty(lx, lx, lx, nelv)
771 real(kind=rp), intent(in) :: wa(lx + 2, lx + 2, lx + 2, nelv)
772 real(kind=rp), intent(in) :: dphidxi(lx, lx)
773
774 call gradient_jump_penalty_finalize_hex(penalty, wa, dphidxi, lx, nelv)
775
776 end subroutine gradient_jump_penalty_finalize
777
785 subroutine gradient_jump_penalty_finalize_hex(penalty, wa, dphidxi, lx, nelv)
786 integer, intent(in) :: lx, nelv
787 real(kind=rp), intent(inout) :: penalty(lx, lx, lx, nelv)
788 real(kind=rp), intent(in) :: wa(lx + 2, lx + 2, lx + 2, nelv)
789 real(kind=rp), intent(in) :: dphidxi(lx, lx)
790
791 integer :: i, j, k
792
793 do i = 1, lx
794 do j = 1, lx
795 do k = 1, lx
796 penalty(i, j, k, :) = &
797 wa(1, j + 1, k + 1, :) * &
798 dphidxi(1, i) + &
799 wa(lx + 2, j + 1, k + 1, :) * &
800 dphidxi(lx, i) + &
801 wa(i + 1, 1, k + 1, :) * &
802 dphidxi(1, j) + &
803 wa(i + 1, lx + 2, k + 1, :) * &
804 dphidxi(lx, j) + &
805 wa(i + 1, j + 1, 1, :) * &
806 dphidxi(1, k) + &
807 wa(i + 1, j + 1, lx + 2, :) * &
808 dphidxi(lx, k)
809 end do
810 end do
811 end do
812
813 end subroutine gradient_jump_penalty_finalize_hex
814
817 subroutine g_compute(this, s)
818 class(gradient_jump_penalty_t), intent(inout) :: this
819 type(field_t), intent(in) :: s
820
821 call dudxyz(this%grad1, s%x, this%coef%drdx, &
822 this%coef%dsdx, this%coef%dtdx, this%coef)
823 call dudxyz(this%grad2, s%x, this%coef%drdy, &
824 this%coef%dsdy, this%coef%dtdy, this%coef)
825 call dudxyz(this%grad3, s%x, this%coef%drdz, &
826 this%coef%dsdz, this%coef%dtdz, this%coef)
827
828 if (neko_bcknd_device .eq. 1) then
829 call device_pick_facet_value_hex(this%flux1_d, this%grad1_d, &
830 this%lx, this%coef%msh%nelv)
831 call device_pick_facet_value_hex(this%flux2_d, this%grad2_d, &
832 this%lx, this%coef%msh%nelv)
833 call device_pick_facet_value_hex(this%flux3_d, this%grad3_d, &
834 this%lx, this%coef%msh%nelv)
835 call device_col2(this%flux1_d, this%n1_d, this%n_large)
836 call device_col2(this%flux2_d, this%n2_d, this%n_large)
837 call device_col2(this%flux3_d, this%n3_d, this%n_large)
838 call device_add3s2(this%G_d, this%flux1_d, this%flux2_d, &
839 1.0_rp, 1.0_rp, this%n_large)
840 call device_add2(this%G_d, this%flux3_d, this%n_large)
841 else
842 call pick_facet_value_hex(this%flux1, this%grad1, &
843 this%lx, this%coef%msh%nelv)
844 call pick_facet_value_hex(this%flux2, this%grad2, &
845 this%lx, this%coef%msh%nelv)
846 call pick_facet_value_hex(this%flux3, this%grad3, &
847 this%lx, this%coef%msh%nelv)
848 call col2(this%flux1, this%n1, this%n_large)
849 call col2(this%flux2, this%n2, this%n_large)
850 call col2(this%flux3, this%n3, this%n_large)
851 call add3(this%G, this%flux1, this%flux2, this%n_large)
852 call add2(this%G, this%flux3, this%n_large)
853 end if
854
855 call this%gs_GJP%op(this%G, this%n_large, gs_op_add)
856
857 end subroutine g_compute
858
863 subroutine absvolflux_compute(this, u, v, w)
864 class(gradient_jump_penalty_t), intent(inout) :: this
865 type(field_t), intent(in) :: u, v, w
866
867 integer :: i
868
869 if (neko_bcknd_device .eq. 1) then
870 call device_pick_facet_value_hex(this%volflux1_d, u%x_d, this%lx, &
871 this%coef%msh%nelv)
872 call device_pick_facet_value_hex(this%volflux2_d, v%x_d, this%lx, &
873 this%coef%msh%nelv)
874 call device_pick_facet_value_hex(this%volflux3_d, w%x_d, this%lx, &
875 this%coef%msh%nelv)
876 call device_col2(this%volflux1_d, this%n1_d, this%n_large)
877 call device_col2(this%volflux2_d, this%n2_d, this%n_large)
878 call device_col2(this%volflux3_d, this%n3_d, this%n_large)
879 call device_add3s2(this%absvolflux_d, this%volflux1_d, &
880 this%volflux2_d, 1.0_rp, 1.0_rp, this%n_large)
881 call device_add2(this%absvolflux_d, this%volflux3_d, this%n_large)
882 call device_absval(this%absvolflux_d, this%n_large)
883 else
884 call pick_facet_value_hex(this%volflux1, u%x, &
885 this%lx, this%coef%msh%nelv)
886 call pick_facet_value_hex(this%volflux2, v%x, &
887 this%lx, this%coef%msh%nelv)
888 call pick_facet_value_hex(this%volflux3, w%x, &
889 this%lx, this%coef%msh%nelv)
890 call col2(this%volflux1, this%n1, this%n_large)
891 call col2(this%volflux2, this%n2, this%n_large)
892 call col2(this%volflux3, this%n3, this%n_large)
893 call add3(this%absvolflux, this%volflux1, this%volflux2, this%n_large)
894 call add2(this%absvolflux, this%volflux3, this%n_large)
895 call absval(this%absvolflux, this%n_large)
896 end if
897
898 end subroutine absvolflux_compute
899
905 subroutine pick_facet_value_hex(f_facet, f_field, lx, nelv)
906 integer, intent(in) :: lx, nelv
907 real(kind=rp), intent(in) :: f_field(lx, lx, lx, nelv)
908 real(kind=rp), intent(inout) :: f_facet(lx + 2, lx + 2, lx + 2, nelv)
909
910 call copy(f_facet(1, 2: lx + 1, 2: lx + 1, :), &
911 f_field(1, :, :, :), lx * lx * nelv)
912 call copy(f_facet(lx + 2, 2: lx + 1, 2: lx + 1, :), &
913 f_field(lx, :, :, :), lx * lx * nelv)
914 call copy(f_facet(2: lx + 1, 1, 2: lx + 1, :), &
915 f_field(:, 1, :, :), lx * lx * nelv)
916 call copy(f_facet(2: lx + 1, lx + 2, 2: lx + 1, :), &
917 f_field(:, lx, :, :), lx * lx * nelv)
918 call copy(f_facet(2: lx + 1, 2: lx + 1, 1, :), &
919 f_field(:, :, 1, :), lx * lx * nelv)
920 call copy(f_facet(2: lx + 1, 2: lx + 1, lx + 2, :), &
921 f_field(:, :, lx, :), lx * lx * nelv)
922
923 end subroutine pick_facet_value_hex
924
925end module gradient_jump_penalty
device::device_map
Map a Fortran array to a device (allocate and associate)
Definition device.F90:57
device::device_memcpy
Copy data between host and device (or device and device)
Definition device.F90:51
json_utils::json_get_or_default
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Definition json_utils.f90:54
utils::neko_error
Definition utils.f90:42
coefs
Coefficients.
Definition coef.f90:34
device_gradient_jump_penalty
Definition device_gradient_jump_penalty.F90:33
device_math
Definition device_math.F90:33
device
Device abstraction, common interface for various accelerators.
Definition device.F90:34
device::host_to_device
integer, parameter, public host_to_device
Definition device.F90:47
dofmap
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
field
Defines a field.
Definition field.f90:34
gather_scatter
Gather-scatter.
Definition gather_scatter.f90:34
gradient_jump_penalty
Implements gradient_jump_penalty_t.
Definition gradient_jump_penalty.f90:35
gradient_jump_penalty::gradient_jump_penalty_finalize_hex
subroutine gradient_jump_penalty_finalize_hex(penalty, wa, dphidxi, lx, nelv)
Finalizinge the gradient jump penalty term for hexahedral elements. <tau * h^2 * absvolflux * G * phi...
Definition gradient_jump_penalty.f90:786
gradient_jump_penalty::gradient_jump_penalty_compute
subroutine gradient_jump_penalty_compute(this, u, v, w, s)
Compute the gradient jump penalty term.
Definition gradient_jump_penalty.f90:718
gradient_jump_penalty::gradient_jump_penalty_finalize
subroutine gradient_jump_penalty_finalize(penalty, wa, dphidxi, lx, nelv)
Interface of finalizing the gradient jump penalty term. <tau * h^2 * absvolflux * G * phij * phik * d...
Definition gradient_jump_penalty.f90:769
gradient_jump_penalty::absvolflux_compute
subroutine absvolflux_compute(this, u, v, w)
Compute the average of the volumetric flux over facets.
Definition gradient_jump_penalty.f90:864
gradient_jump_penalty::gradient_jump_penalty_perform
subroutine gradient_jump_penalty_perform(this, f)
Assign the gradient jump penalty term.
Definition gradient_jump_penalty.f90:748
gradient_jump_penalty::facet_factor_init
subroutine facet_factor_init(this)
Initialize the facet factor array.
Definition gradient_jump_penalty.f90:498
gradient_jump_penalty::g_compute
subroutine g_compute(this, s)
Compute the average of the flux over facets.
Definition gradient_jump_penalty.f90:818
gradient_jump_penalty::eval_h2_hex
subroutine eval_h2_hex(h2_el, n, i, dm, coef)
Evaluate h^2 for each element for hexahedral mesh.
Definition gradient_jump_penalty.f90:335
gradient_jump_penalty::dist2_quadrature_hex
real(kind=rp) function dist2_quadrature_hex(l, k, j, i, n, dm, coef)
Definition gradient_jump_penalty.f90:376
gradient_jump_penalty::gradient_jump_penalty_init
subroutine gradient_jump_penalty_init(this, params, dofmap, coef, a, b)
Constructor.
Definition gradient_jump_penalty.f90:146
gradient_jump_penalty::pick_facet_value_hex
subroutine pick_facet_value_hex(f_facet, f_field, lx, nelv)
Pick facet values of a field.
Definition gradient_jump_penalty.f90:906
gradient_jump_penalty::gradient_jump_penalty_free
subroutine gradient_jump_penalty_free(this)
Destructor for the gradient_jump_penalty_t class.
Definition gradient_jump_penalty.f90:585
gs_ops
Defines Gather-scatter operations.
Definition gs_ops.f90:34
gs_ops::gs_op_add
integer, parameter, public gs_op_add
Definition gs_ops.f90:36
hex
Defines a hexahedron element.
Definition hex.f90:34
json_utils
Utilities for retrieving parameters from the case files.
Definition json_utils.f90:34
math
Definition math.f90:60
neko_config
Build configurations.
Definition neko_config.f90:34
neko_config::neko_bcknd_device
integer, parameter neko_bcknd_device
Definition neko_config.f90:44
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
operators
Operators.
Definition operators.f90:34
operators::dudxyz
subroutine, public dudxyz(du, u, dr, ds, dt, coef)
Compute derivative of a scalar field along a single direction.
Definition operators.f90:76
point
Implements a point.
Definition point.f90:35
quad
Defines a quadrilateral element.
Definition quad.f90:34
space
Defines a function space.
Definition space.f90:34
space::gll
integer, parameter, public gll
Definition space.f90:48
utils
Utilities.
Definition utils.f90:35
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
dofmap::dofmap_t
Definition dofmap.f90:53
element::element_t
Base type for an element.
Definition element.f90:44
field::field_t
Definition field.f90:46
gather_scatter::gs_t
Definition gather_scatter.f90:58
gradient_jump_penalty::gradient_jump_penalty_t
Implements the gradient jump penalty.
Definition gradient_jump_penalty.f90:62
hex::hex_t
Hexahedron element.
Definition hex.f90:63
point::point_t
A point in with coordinates .
Definition point.f90:43
quad::quad_t
Quadrilateral element.
Definition quad.f90:58
space::space_t
The function space for the SEM solution fields.
Definition space.f90:62