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