Neko 1.99.2
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!
36 use num_types, only : rp
37 use utils, only : neko_error
39 use json_module, only : json_file
40 use math, only: add2, col2, col3, invcol2, add3, copy, absval
41 use point, only : point_t
42 use field, only : field_t
43 use dofmap , only : dofmap_t
45 use coefs, only : coef_t
46 use element, only : element_t
47 use hex, only : hex_t
48 use quad, only : quad_t
49 use gs_ops, only : gs_op_add
50 use space, only : space_t, gll
51 use gather_scatter, only : gs_t
57 use source_term, only : source_term_t
58 use field_list, only : field_list_t
59 use registry, only : neko_registry
60 use time_state, only : time_state_t
61 use operators, only : dudxyz
62 use, intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr, c_associated
63
64 implicit none
65 private
66
69 type, public, extends(source_term_t):: gradient_jump_penalty_t
71 type(field_list_t) :: s_fields
73 type(field_t), pointer :: u
75 type(field_t), pointer :: v
77 type(field_t), pointer :: w
78
80 real(kind=rp) :: tau
82 integer :: p
84 integer :: lx
86 real(kind=rp), allocatable, dimension(:, :, :, :) :: penalty
87 type(c_ptr) :: penalty_d = c_null_ptr
89 real(kind=rp), allocatable, dimension(:, :, :, :) :: penalty_facet
90 type(c_ptr) :: penalty_facet_d = c_null_ptr
92 real(kind=rp), allocatable, dimension(:, :, :, :) :: grad1, &
93 grad2, grad3
94 type(c_ptr) :: grad1_d = c_null_ptr
95 type(c_ptr) :: grad2_d = c_null_ptr
96 type(c_ptr) :: grad3_d = c_null_ptr
98 real(kind=rp), allocatable, dimension(:, :, :, :) :: g
99 type(c_ptr) :: g_d = c_null_ptr
101 real(kind=rp), allocatable, dimension(:, :, :, :) :: flux1, flux2, flux3
102 type(c_ptr) :: flux1_d = c_null_ptr
103 type(c_ptr) :: flux2_d = c_null_ptr
104 type(c_ptr) :: flux3_d = c_null_ptr
106 real(kind=rp), allocatable, dimension(:, :, :, :) :: volflux1, &
107 volflux2, volflux3
108 type(c_ptr) :: volflux1_d = c_null_ptr
109 type(c_ptr) :: volflux2_d = c_null_ptr
110 type(c_ptr) :: volflux3_d = c_null_ptr
112 real(kind=rp), allocatable, dimension(:, :, :, :) :: absvolflux
113 type(c_ptr) :: absvolflux_d = c_null_ptr
115 real(kind=rp), allocatable, dimension(:, :, :, :) :: n1, n2, n3
116 type(c_ptr) :: n1_d = c_null_ptr
117 type(c_ptr) :: n2_d = c_null_ptr
118 type(c_ptr) :: n3_d = c_null_ptr
120 real(kind=rp), allocatable, dimension(:, :, :, :) :: facet_factor
121 type(c_ptr) :: facet_factor_d = c_null_ptr
123 integer, allocatable :: n_facet(:)
124 integer :: n_facet_max
126 real(kind=rp), allocatable, dimension(:, :, :, :) :: h2
128 real(kind=rp), allocatable :: dphidxi(:, :)
129 type(c_ptr) :: dphidxi_d = c_null_ptr
131 type(space_t) :: xh_gjp
132 type(dofmap_t) :: dm_gjp
133 type(gs_t) :: gs_gjp
135 integer :: n
137 integer :: n_large
138
139 contains
141 procedure, pass(this) :: init => gradient_jump_penalty_init
143 procedure, pass(this) :: init_from_components => &
146 procedure, pass(this) :: free => gradient_jump_penalty_free
148 procedure, pass(this) :: compute_single => &
151 procedure, pass(this) :: compute_ => gradient_jump_penalty_compute
152
154
155contains
162 subroutine gradient_jump_penalty_init(this, json, fields, coef, variable_name)
163 implicit none
164 class(gradient_jump_penalty_t), intent(inout) :: this
165 type(json_file), intent(inout) :: json
166 type(coef_t), target, intent(in) :: coef
167 character(len=*), intent(in) :: variable_name
168 type(field_list_t), intent(in), target :: fields
169
170 real(kind=rp) :: start_time, end_time
171 real(kind=rp) :: a, b
172
173 call json_get_or_default(json, "start_time", start_time, 0.0_rp)
174 call json_get_or_default(json, "end_time", end_time, huge(0.0_rp))
175
176 if ((coef%Xh%lx - 1) .eq. 1) then
177 call json_get_or_default(json, 'tau', a, 0.02_rp)
178 b = 0.0_rp
179 else
180 call json_get_or_default(json, 'scaling_factor', a, 0.8_rp)
181 call json_get_or_default(json, 'scaling_exponent', b, 4.0_rp)
182 end if
183
184 call this%init_from_components(coef, fields, start_time, end_time, a, b, &
185 variable_name)
186
187 end subroutine gradient_jump_penalty_init
188
196 subroutine gradient_jump_penalty_init_from_components(this, coef, fields, &
197 start_time, end_time, a, b, variable_name)
198 class(gradient_jump_penalty_t), intent(inout) :: this
199 type(coef_t), target, intent(in) :: coef
200 type(field_list_t), intent(in), target :: fields
201 real(kind=rp), intent(in) :: start_time
202 real(kind=rp), intent(in) :: end_time
203 real(kind=rp), intent(in) :: a, b
204 character(len=*), intent(in) :: variable_name
205
206 integer :: i, j, k, l
207 real(kind=rp), allocatable :: zg(:) ! Quadrature points
208 real(kind=rp) :: normal(3)
209
210 call this%free()
211
212 call this%init_base(fields, coef, start_time, end_time)
213
214 this%u => neko_registry%get_field("u")
215 this%v => neko_registry%get_field("v")
216 this%w => neko_registry%get_field("w")
217
218 if (fields%size() .eq. 1) then
219 call this%s_fields%init(1)
220 call this%s_fields%assign(1, &
221 neko_registry%get_field(variable_name))
222 else if (fields%size() .eq. 3) then
223 call this%s_fields%init(3)
224 call this%s_fields%assign(1, this%u)
225 call this%s_fields%assign(2, this%v)
226 call this%s_fields%assign(3, this%w)
227 else
228 call neko_error("The GJP source assumes either 3 or 1 RHS fields.")
229 end if
230
231
232
233
234 this%p = coef%dof%Xh%lx - 1
235 this%lx = coef%dof%Xh%lx
236
237 if (this%p .gt. 1) then
238 this%tau = -a * (this%p + 1) ** (-b)
239 else
240 this%tau = -a
241 end if
242
243 this%n = this%lx ** 3 * this%coef%msh%nelv
244 this%n_large = (this%lx + 2) ** 3 * this%coef%msh%nelv
245
246 allocate(this%n_facet(this%coef%msh%nelv))
247 do i = 1, this%coef%msh%nelv
248 select type (ep => this%coef%msh%elements(i)%e)
249 type is (hex_t)
250 this%n_facet(i) = 6
251 type is (quad_t)
252 call neko_error("Only Hexahedral element is &
253 &supported now for gradient jump penalty")
254 end select
255 end do
256 this%n_facet_max = maxval(this%n_facet)
257
258 allocate(this%h2(this%lx + 2, this%lx + 2, &
259 this%lx + 2, this%coef%msh%nelv))
260
261 do i = 1, this%coef%msh%nelv
262 select type (ep => this%coef%msh%elements(i)%e)
263 type is (hex_t)
264 call eval_h2_hex(this%h2(:, :, :, i), this%lx, i, this%coef)
265 type is (quad_t)
266 call neko_error("Gradient jump penalty error: mesh size &
267 &evaluation is not supported for quad_t")
268 end select
269 end do
270
271 allocate(zg(this%lx))
272 allocate(this%dphidxi(this%lx, this%lx))
273
274 zg = coef%Xh%zg(:,1)
275 do i = 1, coef%Xh%lx
276 do j = 1, coef%Xh%lx
277 this%dphidxi(j,i) = this%coef%Xh%dx(j,i)
278 end do
279 end do
280
281 allocate(this%penalty(this%lx, this%lx, this%lx, this%coef%msh%nelv))
282 allocate(this%grad1(this%lx, this%lx, this%lx, this%coef%msh%nelv))
283 allocate(this%grad2(this%lx, this%lx, this%lx, this%coef%msh%nelv))
284 allocate(this%grad3(this%lx, this%lx, this%lx, this%coef%msh%nelv))
285
286 allocate(this%penalty_facet(this%lx + 2, this%lx + 2, &
287 this%lx + 2, this%coef%msh%nelv))
288 allocate(this%G(this%lx + 2, this%lx + 2, &
289 this%lx + 2, this%coef%msh%nelv))
290 allocate(this%flux1(this%lx + 2, this%lx + 2, &
291 this%lx + 2, this%coef%msh%nelv))
292 allocate(this%flux2(this%lx + 2, this%lx + 2, &
293 this%lx + 2, this%coef%msh%nelv))
294 allocate(this%flux3(this%lx + 2, this%lx + 2, &
295 this%lx + 2, this%coef%msh%nelv))
296 allocate(this%volflux1(this%lx + 2, this%lx + 2, &
297 this%lx + 2, this%coef%msh%nelv))
298 allocate(this%volflux2(this%lx + 2, this%lx + 2, &
299 this%lx + 2, this%coef%msh%nelv))
300 allocate(this%volflux3(this%lx + 2, this%lx + 2, &
301 this%lx + 2, this%coef%msh%nelv))
302 allocate(this%absvolflux(this%lx + 2, this%lx + 2, &
303 this%lx + 2, this%coef%msh%nelv))
304 allocate(this%n1(this%lx + 2, this%lx + 2, &
305 this%lx + 2, this%coef%msh%nelv))
306 allocate(this%n2(this%lx + 2, this%lx + 2, &
307 this%lx + 2, this%coef%msh%nelv))
308 allocate(this%n3(this%lx + 2, this%lx + 2, &
309 this%lx + 2, this%coef%msh%nelv))
310
311 ! Extract facets' normals
312 do i = 1, this%coef%msh%nelv
313 do j = 1, 6 ! for hexahedral elements
314 do k = 1, this%lx
315 do l = 1, this%lx
316 select case (j)
317 case (1)
318 normal = this%coef%get_normal(1, l, k, i, j)
319 this%n1(1, l + 1, k + 1, i) = normal(1)
320 this%n2(1, l + 1, k + 1, i) = normal(2)
321 this%n3(1, l + 1, k + 1, i) = normal(3)
322 case (2)
323 normal = this%coef%get_normal(1, l, k, i, j)
324 this%n1(this%lx + 2, l + 1, k + 1, i) = normal(1)
325 this%n2(this%lx + 2, l + 1, k + 1, i) = normal(2)
326 this%n3(this%lx + 2, l + 1, k + 1, i) = normal(3)
327 case (3)
328 normal = this%coef%get_normal(l, 1, k, i, j)
329 this%n1(l + 1, 1, k + 1, i) = normal(1)
330 this%n2(l + 1, 1, k + 1, i) = normal(2)
331 this%n3(l + 1, 1, k + 1, i) = normal(3)
332 case (4)
333 normal = this%coef%get_normal(l, 1, k, i, j)
334 this%n1(l + 1, this%lx + 2, k + 1, i) = normal(1)
335 this%n2(l + 1, this%lx + 2, k + 1, i) = normal(2)
336 this%n3(l + 1, this%lx + 2, k + 1, i) = normal(3)
337 case (5)
338 normal = this%coef%get_normal(l, k, 1, i, j)
339 this%n1(l + 1, k + 1, 1, i) = normal(1)
340 this%n2(l + 1, k + 1, 1, i) = normal(2)
341 this%n3(l + 1, k + 1, 1, i) = normal(3)
342 case (6)
343 normal = this%coef%get_normal(l, k, 1, i, j)
344 this%n1(l + 1, k + 1, this%lx + 2, i) = normal(1)
345 this%n2(l + 1, k + 1, this%lx + 2, i) = normal(2)
346 this%n3(l + 1, k + 1, this%lx + 2, i) = normal(3)
347 case default
348 call neko_error("The face index is not correct")
349 end select
350 end do
351 end do
352 end do
353 end do
354
355 ! Assemble facet factor
356 call facet_factor_init(this)
357
358 ! Initialize Gather-Scatter
359 call this%Xh_GJP%init(gll, this%lx+2, this%lx+2, this%lx+2)
360 call this%dm_GJP%init(this%coef%msh, this%Xh_GJP)
361 call this%gs_GJP%init(this%dm_GJP)
362
363 ! Initialize pointers for device
364 if (neko_bcknd_device .eq. 1) then
365 call device_map(this%dphidxi, this%dphidxi_d, &
366 this%lx * this%lx)
367 call device_map(this%penalty, this%penalty_d, this%n)
368 call device_map(this%grad1, this%grad1_d, this%n)
369 call device_map(this%grad2, this%grad2_d, this%n)
370 call device_map(this%grad3, this%grad3_d, this%n)
371
372 call device_map(this%penalty_facet, this%penalty_facet_d, this%n_large)
373 call device_map(this%G, this%G_d, this%n_large)
374 call device_map(this%flux1, this%flux1_d, this%n_large)
375 call device_map(this%flux2, this%flux2_d, this%n_large)
376 call device_map(this%flux3, this%flux3_d, this%n_large)
377
378 call device_map(this%volflux1, this%volflux1_d, this%n_large)
379 call device_map(this%volflux2, this%volflux2_d, this%n_large)
380 call device_map(this%volflux3, this%volflux3_d, this%n_large)
381 call device_map(this%absvolflux, this%absvolflux_d, this%n_large)
382
383 call device_map(this%n1, this%n1_d, this%n_large)
384 call device_map(this%n2, this%n2_d, this%n_large)
385 call device_map(this%n3, this%n3_d, this%n_large)
386 call device_map(this%facet_factor, this%facet_factor_d, this%n_large)
387
388 call device_memcpy(this%dphidxi, this%dphidxi_d, &
389 this%lx * this%lx, &
390 host_to_device, sync = .false.)
391 call device_memcpy(this%n1, this%n1_d, this%n_large, &
392 host_to_device, sync = .false.)
393 call device_memcpy(this%n2, this%n2_d, this%n_large, &
394 host_to_device, sync = .false.)
395 call device_memcpy(this%n3, this%n3_d, this%n_large, &
396 host_to_device, sync = .false.)
397 call device_memcpy(this%facet_factor, this%facet_factor_d, this%n_large,&
398 host_to_device, sync = .false.)
399
400 end if
401
403
407 subroutine eval_h2_hex(h2_el, n, i, coef)
408 integer, intent(in) :: n, i
409 type(coef_t), pointer, intent(in) :: coef
410 real(kind=rp), intent(inout) :: h2_el(n + 2, n + 2, n + 2)
411
412 type(dofmap_t), pointer :: dm
413 integer :: j, k, l
414
415 dm => coef%dof
416
417 h2_el = 0.0_rp
418
419 do j = 1, 6
420 do k = 1, n
421 do l = 1, n
422 select case (j)
423 case (1)
424 h2_el(1, l + 1, k + 1) = &
425 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
426 case (2)
427 h2_el(n + 2, l + 1, k + 1) = &
428 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
429 case (3)
430 h2_el(l + 1, 1, k + 1) = &
431 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
432 case (4)
433 h2_el(l + 1, n + 2, k + 1) = &
434 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
435 case (5)
436 h2_el(l + 1, k + 1, 1) = &
437 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
438 case (6)
439 h2_el(l + 1, k + 1, n + 2) = &
440 dist2_quadrature_hex(l, k, j, i, n, dm, coef)
441 case default
442 call neko_error("The face index is not correct")
443 end select
444 end do
445 end do
446 end do
447
448 end subroutine eval_h2_hex
449
450 function dist2_quadrature_hex(l, k, j, i, n, dm, coef) result(dist2)
451 integer, intent(in) :: l, k, j, i, n
452 type(dofmap_t), pointer, intent(in) :: dm
453 type(coef_t), pointer, intent(in) :: coef
454 real(kind=rp) :: dist2, dist_1, dist_2
455
456 real(kind=rp) :: x1, y1, z1, x2, y2, z2
457 real(kind=rp) :: normal1(3), normal2(3)
458 real(kind=rp) :: n11, n12, n13, n21, n22, n23
459 real(kind=rp) :: v1, v2, v3
460
461 dist2 = 0.0_rp
462 select case (j)
463 case (1)
464 normal1 = coef%get_normal(1, l, k, i, 1)
465 n11 = normal1(1)
466 n12 = normal1(2)
467 n13 = normal1(3)
468 x1 = dm%x(1, l, k, i)
469 y1 = dm%y(1, l, k, i)
470 z1 = dm%z(1, l, k, i)
471 normal2 = coef%get_normal(1, l, k, i, 2)
472 n21 = normal2(1)
473 n22 = normal2(2)
474 n23 = normal2(3)
475 x2 = dm%x(n, l, k, i)
476 y2 = dm%y(n, l, k, i)
477 z2 = dm%z(n, l, k, i)
478 case (2)
479 ! now the facet pair share the same value for h
480 ! but just let it be here for furture possible changes
481 normal1 = coef%get_normal(1, l, k, i, 2)
482 n11 = normal1(1)
483 n12 = normal1(2)
484 n13 = normal1(3)
485 x1 = dm%x(n, l, k, i)
486 y1 = dm%y(n, l, k, i)
487 z1 = dm%z(n, l, k, i)
488 normal2 = coef%get_normal(1, l, k, i, 1)
489 n21 = normal2(1)
490 n22 = normal2(2)
491 n23 = normal2(3)
492 x2 = dm%x(1, l, k, i)
493 y2 = dm%y(1, l, k, i)
494 z2 = dm%z(1, l, k, i)
495 case (3)
496 normal1 = coef%get_normal(1, l, k, i, 3)
497 n11 = normal1(1)
498 n12 = normal1(2)
499 n13 = normal1(3)
500 x1 = dm%x(l, 1, k, i)
501 y1 = dm%y(l, 1, k, i)
502 z1 = dm%z(l, 1, k, i)
503 normal2 = coef%get_normal(1, l, k, i, 4)
504 n21 = normal2(1)
505 n22 = normal2(2)
506 n23 = normal2(3)
507 x2 = dm%x(l, n, k, i)
508 y2 = dm%y(l, n, k, i)
509 z2 = dm%z(l, n, k, i)
510 case (4)
511 normal1 = coef%get_normal(1, l, k, i, 4)
512 n11 = normal1(1)
513 n12 = normal1(2)
514 n13 = normal1(3)
515 x1 = dm%x(l, n, k, i)
516 y1 = dm%y(l, n, k, i)
517 z1 = dm%z(l, n, k, i)
518 normal2 = coef%get_normal(1, l, k, i, 3)
519 n21 = normal2(1)
520 n22 = normal2(2)
521 n23 = normal2(3)
522 x2 = dm%x(l, 1, k, i)
523 y2 = dm%y(l, 1, k, i)
524 z2 = dm%z(l, 1, k, i)
525 case (5)
526 normal1 = coef%get_normal(1, l, k, i, 5)
527 n11 = normal1(1)
528 n12 = normal1(2)
529 n13 = normal1(3)
530 x1 = dm%x(l, k, 1, i)
531 y1 = dm%y(l, k, 1, i)
532 z1 = dm%z(l, k, 1, i)
533 normal2 = coef%get_normal(1, l, k, i, 6)
534 n21 = normal2(1)
535 n22 = normal2(2)
536 n23 = normal2(3)
537 x2 = dm%x(l, k, n, i)
538 y2 = dm%y(l, k, n, i)
539 z2 = dm%z(l, k, n, i)
540 case (6)
541 normal1 = coef%get_normal(1, l, k, i, 6)
542 n11 = normal1(1)
543 n12 = normal1(2)
544 n13 = normal1(3)
545 x1 = dm%x(l, k, n, i)
546 y1 = dm%y(l, k, n, i)
547 z1 = dm%z(l, k, n, i)
548 normal2 = coef%get_normal(1, l, k, i, 5)
549 n21 = normal2(1)
550 n22 = normal2(2)
551 n23 = normal2(3)
552 x2 = dm%x(l, k, 1, i)
553 y2 = dm%y(l, k, 1, i)
554 z2 = dm%z(l, k, 1, i)
555 case default
556 call neko_error("The face index is not correct")
557 end select
558
559 ! get the vector from the quadrature point to the one on the other side
560 v1 = x2 - x1
561 v2 = y2 - y1
562 v3 = z2 - z1
563 ! Project onto tabsvolflhe facet-normal direction of the point
564 dist_1 = v1*n11 + v2*n12 + v3*n13
565 dist_2 = - (v1*n21 + v2*n22 + v3*n23)
566
567 dist2 = ((dist_1 + dist_2)/2.0_rp)*((dist_1 + dist_2)/2.0_rp)
568
569 end function dist2_quadrature_hex
570
572 subroutine facet_factor_init(this)
573 class(gradient_jump_penalty_t), intent(inout) :: this
574 integer :: i, j, k, l
575 real(kind=rp) :: area_tmp
576
577 allocate(this%facet_factor(this%lx + 2, this%lx + 2, &
578 this%lx + 2, this%coef%msh%nelv))
579
580 associate(facet_factor => this%facet_factor, &
581 coef => this%coef, &
582 lx => this%lx, &
583 nelv => this%coef%msh%nelv, &
584 jacinv => this%coef%jacinv, h2 => this%h2, &
585 tau => this%tau, n1 => this%n1, &
586 n2 => this%n2, n3 => this%n3)
587
588 do i = 1, nelv
589 do j = 1, 6 ! for hexahedral elementsh2
590 do k = 1, lx
591 do l = 1, lx
592 select case (j)
593 case (1)
594 area_tmp = coef%get_area(1, l, k, i, j)
595 facet_factor(1, l + 1, k + 1, i) = area_tmp * tau * &
596 h2(1, l + 1, k + 1, i) * &
597 (n1(1, l + 1, k + 1, i) * coef%drdx(1, l, k, i) + &
598 n2(1, l + 1, k + 1, i) * coef%drdy(1, l, k, i) + &
599 n3(1, l + 1, k + 1, i) * coef%drdz(1, l, k, i) ) &
600 * jacinv(1, l, k, i)
601 case (2)
602 area_tmp = coef%get_area(1, l, k, i, j)
603 facet_factor(lx + 2, l + 1, k + 1, i) = area_tmp * tau * &
604 h2(lx + 2, l + 1, k + 1, i) * &
605 (n1(lx + 2, l + 1, k + 1, i) * &
606 coef%drdx(lx, l, k, i) + &
607 n2(lx + 2, l + 1, k + 1, i) * &
608 coef%drdy(lx, l, k, i) + &
609 n3(lx + 2, l + 1, k + 1, i) * &
610 coef%drdz(lx, l, k, i) ) &
611 * jacinv(lx, l, k, i)
612 case (3)
613 area_tmp = coef%get_area(l, 1, k, i, j)
614 facet_factor(l + 1, 1, k + 1, i) = area_tmp * tau * &
615 h2(l + 1, 1, k + 1, i) * &
616 (n1(l + 1, 1, k + 1, i) * coef%dsdx(l, 1, k, i) + &
617 n2(l + 1, 1, k + 1, i) * coef%dsdy(l, 1, k, i) + &
618 n3(l + 1, 1, k + 1, i) * coef%dsdz(l, 1, k, i) ) &
619 * jacinv(l, 1, k, i)
620 case (4)
621 area_tmp = coef%get_area(l, 1, k, i, j)
622 facet_factor(l + 1, lx + 2, k + 1, i) = area_tmp * tau * &
623 h2(l + 1, lx + 2, k + 1, i) * &
624 (n1(l + 1, lx + 2, k + 1, i) * &
625 coef%dsdx(l, lx, k, i) + &
626 n2(l + 1, lx + 2, k + 1, i) * &
627 coef%dsdy(l, lx, k, i) + &
628 n3(l + 1, lx + 2, k + 1, i) * &
629 coef%dsdz(l, lx, k, i) ) &
630 * jacinv(l, lx, k, i)
631 case (5)
632 area_tmp = coef%get_area(l, k, 1, i, j)
633 facet_factor(l + 1, k + 1, 1, i) = area_tmp * tau * &
634 h2(l + 1, k + 1, 1, i) * &
635 (n1(l + 1, k + 1, 1, i) * coef%dtdx(l, k, 1, i) + &
636 n2(l + 1, k + 1, 1, i) * coef%dtdy(l, k, 1, i) + &
637 n3(l + 1, k + 1, 1, i) * coef%dtdz(l, k, 1, i) ) &
638 * jacinv(l, k, 1, i)
639 case (6)
640 area_tmp = coef%get_area(l, k, 1, i, j)
641 facet_factor(l + 1, k + 1, lx + 2, i) = area_tmp * tau * &
642 h2(l + 1, k + 1, lx + 2, i) * &
643 ( &
644 n1(l + 1, k + 1, lx + 2, i) * &
645 coef%dtdx(l, k, lx, i) + &
646 n2(l + 1, k + 1, lx + 2, i) * &
647 coef%dtdy(l, k, lx, i) + &
648 n3(l + 1, k + 1, lx + 2, i) * &
649 coef%dtdz(l, k, lx, i) &
650 ) &
651 * jacinv(l, k, lx, i)
652 case default
653 call neko_error("The face index is not correct")
654 end select
655 end do
656 end do
657 end do
658 end do
659
660 end associate
661 end subroutine facet_factor_init
662
665 implicit none
666 class(gradient_jump_penalty_t), intent(inout) :: this
667
668 call this%free_base
669
670 if (c_associated(this%dphidxi_d)) then
671 call device_free(this%dphidxi_d)
672 end if
673 if (c_associated(this%penalty_d)) then
674 call device_free(this%penalty_d)
675 end if
676 if (c_associated(this%grad1_d)) then
677 call device_free(this%grad1_d)
678 end if
679 if (c_associated(this%grad2_d)) then
680 call device_free(this%grad2_d)
681 end if
682 if (c_associated(this%grad3_d)) then
683 call device_free(this%grad3_d)
684 end if
685 if (c_associated(this%penalty_facet_d)) then
686 call device_free(this%penalty_facet_d)
687 end if
688 if (c_associated(this%G_d)) then
689 call device_free(this%G_d)
690 end if
691 if (c_associated(this%flux1_d)) then
692 call device_free(this%flux1_d)
693 end if
694 if (c_associated(this%flux2_d)) then
695 call device_free(this%flux2_d)
696 end if
697 if (c_associated(this%flux3_d)) then
698 call device_free(this%flux3_d)
699 end if
700 if (c_associated(this%volflux1_d)) then
701 call device_free(this%volflux1_d)
702 end if
703 if (c_associated(this%volflux2_d)) then
704 call device_free(this%volflux2_d)
705 end if
706 if (c_associated(this%volflux3_d)) then
707 call device_free(this%volflux3_d)
708 end if
709 if (c_associated(this%absvolflux_d)) then
710 call device_free(this%absvolflux_d)
711 end if
712 if (c_associated(this%n1_d)) then
713 call device_free(this%n1_d)
714 end if
715 if (c_associated(this%n2_d)) then
716 call device_free(this%n2_d)
717 end if
718 if (c_associated(this%n3_d)) then
719 call device_free(this%n3_d)
720 end if
721 if (c_associated(this%facet_factor_d)) then
722 call device_free(this%facet_factor_d)
723 end if
724
725 if (allocated(this%penalty)) then
726 deallocate(this%penalty)
727 end if
728 if (allocated(this%grad1)) then
729 deallocate(this%grad1)
730 end if
731 if (allocated(this%grad2)) then
732 deallocate(this%grad2)
733 end if
734 if (allocated(this%grad3)) then
735 deallocate(this%grad3)
736 end if
737 if (allocated(this%h2)) then
738 deallocate(this%h2)
739 end if
740 if (allocated(this%n_facet)) then
741 deallocate(this%n_facet)
742 end if
743 if (allocated(this%dphidxi)) then
744 deallocate(this%dphidxi)
745 end if
746 if (allocated(this%penalty_facet)) then
747 deallocate(this%penalty_facet)
748 end if
749 if (allocated(this%G)) then
750 deallocate(this%G)
751 end if
752 if (allocated(this%flux1)) then
753 deallocate(this%flux1)
754 end if
755 if (allocated(this%flux2)) then
756 deallocate(this%flux2)
757 end if
758 if (allocated(this%flux3)) then
759 deallocate(this%flux3)
760 end if
761 if (allocated(this%volflux1)) then
762 deallocate(this%volflux1)
763 end if
764 if (allocated(this%volflux2)) then
765 deallocate(this%volflux2)
766 end if
767 if (allocated(this%volflux3)) then
768 deallocate(this%volflux3)
769 end if
770 if (allocated(this%absvolflux)) then
771 deallocate(this%absvolflux)
772 end if
773 if (allocated(this%n1)) then
774 deallocate(this%n1)
775 end if
776 if (allocated(this%n2)) then
777 deallocate(this%n2)
778 end if
779 if (allocated(this%n3)) then
780 deallocate(this%n3)
781 end if
782 if (allocated(this%facet_factor)) then
783 deallocate(this%facet_factor)
784 end if
785
786 nullify(this%u)
787 nullify(this%v)
788 nullify(this%w)
789
790 call this%s_fields%free()
791
792 call this%Xh_GJP%free()
793 call this%gs_GJP%free()
794 call this%dm_GJP%free()
795
796 end subroutine gradient_jump_penalty_free
797
801 class(gradient_jump_penalty_t), intent(inout) :: this
802 type(field_t), intent(in) :: s
803
804 class(element_t), pointer :: ep
805 integer :: i
806
807 call g_compute(this, s)
808 call absvolflux_compute(this, this%u, this%v, this%w)
809
810 if (neko_bcknd_device .eq. 1) then
811 call device_col3(this%penalty_facet_d, this%absvolflux_d, this%G_d, &
812 this%n_large)
813 call device_col2(this%penalty_facet_d, this%facet_factor_d, this%n_large)
814 call device_gradient_jump_penalty_finalize(this%penalty_d, &
815 this%penalty_facet_d, &
816 this%dphidxi_d, &
817 this%lx, this%coef%msh%nelv)
818 else
819 call col3(this%penalty_facet, this%absvolflux, this%G, this%n_large)
820 call col2(this%penalty_facet, this%facet_factor, this%n_large)
821 call gradient_jump_penalty_finalize(this%penalty, this%penalty_facet, &
822 this%dphidxi, &
823 this%lx, this%coef%msh%nelv)
824 end if
825
827
830 subroutine gradient_jump_penalty_compute(this, time)
831 class(gradient_jump_penalty_t), intent(inout) :: this
832 type(time_state_t), intent(in) :: time
833 integer :: i, n_fields, n
834
835 n_fields = this%fields%size()
836 n = this%coef%dof%size()
837
838
839
840 do i = 1, n_fields
841
842 call this%compute_single(this%s_fields%items(i)%ptr)
843
844 if (neko_bcknd_device .eq. 1) then
845 call device_invcol2(this%penalty_d, this%coef%B_d, n)
846 call device_add2(this%fields%x_d(i), this%penalty_d, n)
847 else
848 call invcol2(this%penalty, this%coef%B, n)
849 call add2(this%fields%items(i)%ptr%x, this%penalty, n)
850 end if
851 end do
852
853 end subroutine gradient_jump_penalty_compute
854
862 subroutine gradient_jump_penalty_finalize(penalty, wa, dphidxi, lx, nelv)
863 integer, intent(in) :: lx, nelv
864 real(kind=rp), intent(inout) :: penalty(lx, lx, lx, nelv)
865 real(kind=rp), intent(in) :: wa(lx + 2, lx + 2, lx + 2, nelv)
866 real(kind=rp), intent(in) :: dphidxi(lx, lx)
867
868 call gradient_jump_penalty_finalize_hex(penalty, wa, dphidxi, lx, nelv)
869
870 end subroutine gradient_jump_penalty_finalize
871
879 subroutine gradient_jump_penalty_finalize_hex(penalty, wa, dphidxi, lx, nelv)
880 integer, intent(in) :: lx, nelv
881 real(kind=rp), intent(inout) :: penalty(lx, lx, lx, nelv)
882 real(kind=rp), intent(in) :: wa(lx + 2, lx + 2, lx + 2, nelv)
883 real(kind=rp), intent(in) :: dphidxi(lx, lx)
884
885 integer :: i, j, k
886
887 do i = 1, lx
888 do j = 1, lx
889 do k = 1, lx
890 penalty(i, j, k, :) = &
891 wa(1, j + 1, k + 1, :) * &
892 dphidxi(1, i) + &
893 wa(lx + 2, j + 1, k + 1, :) * &
894 dphidxi(lx, i) + &
895 wa(i + 1, 1, k + 1, :) * &
896 dphidxi(1, j) + &
897 wa(i + 1, lx + 2, k + 1, :) * &
898 dphidxi(lx, j) + &
899 wa(i + 1, j + 1, 1, :) * &
900 dphidxi(1, k) + &
901 wa(i + 1, j + 1, lx + 2, :) * &
902 dphidxi(lx, k)
903 end do
904 end do
905 end do
906
908
911 subroutine g_compute(this, s)
912 class(gradient_jump_penalty_t), intent(inout) :: this
913 type(field_t), intent(in) :: s
914
915 call dudxyz(this%grad1, s%x, this%coef%drdx, &
916 this%coef%dsdx, this%coef%dtdx, this%coef)
917 call dudxyz(this%grad2, s%x, this%coef%drdy, &
918 this%coef%dsdy, this%coef%dtdy, this%coef)
919 call dudxyz(this%grad3, s%x, this%coef%drdz, &
920 this%coef%dsdz, this%coef%dtdz, this%coef)
921
922 if (neko_bcknd_device .eq. 1) then
923 call device_pick_facet_value_hex(this%flux1_d, this%grad1_d, &
924 this%lx, this%coef%msh%nelv)
925 call device_pick_facet_value_hex(this%flux2_d, this%grad2_d, &
926 this%lx, this%coef%msh%nelv)
927 call device_pick_facet_value_hex(this%flux3_d, this%grad3_d, &
928 this%lx, this%coef%msh%nelv)
929 call device_col2(this%flux1_d, this%n1_d, this%n_large)
930 call device_col2(this%flux2_d, this%n2_d, this%n_large)
931 call device_col2(this%flux3_d, this%n3_d, this%n_large)
932 call device_add3s2(this%G_d, this%flux1_d, this%flux2_d, &
933 1.0_rp, 1.0_rp, this%n_large)
934 call device_add2(this%G_d, this%flux3_d, this%n_large)
935 else
936 call pick_facet_value_hex(this%flux1, this%grad1, &
937 this%lx, this%coef%msh%nelv)
938 call pick_facet_value_hex(this%flux2, this%grad2, &
939 this%lx, this%coef%msh%nelv)
940 call pick_facet_value_hex(this%flux3, this%grad3, &
941 this%lx, this%coef%msh%nelv)
942 call col2(this%flux1, this%n1, this%n_large)
943 call col2(this%flux2, this%n2, this%n_large)
944 call col2(this%flux3, this%n3, this%n_large)
945 call add3(this%G, this%flux1, this%flux2, this%n_large)
946 call add2(this%G, this%flux3, this%n_large)
947 end if
948
949 call this%gs_GJP%op(this%G, this%n_large, gs_op_add)
950
951 end subroutine g_compute
952
957 subroutine absvolflux_compute(this, u, v, w)
958 class(gradient_jump_penalty_t), intent(inout) :: this
959 type(field_t), intent(in) :: u, v, w
960
961 integer :: i
962
963 if (neko_bcknd_device .eq. 1) then
964 call device_pick_facet_value_hex(this%volflux1_d, u%x_d, this%lx, &
965 this%coef%msh%nelv)
966 call device_pick_facet_value_hex(this%volflux2_d, v%x_d, this%lx, &
967 this%coef%msh%nelv)
968 call device_pick_facet_value_hex(this%volflux3_d, w%x_d, this%lx, &
969 this%coef%msh%nelv)
970 call device_col2(this%volflux1_d, this%n1_d, this%n_large)
971 call device_col2(this%volflux2_d, this%n2_d, this%n_large)
972 call device_col2(this%volflux3_d, this%n3_d, this%n_large)
973 call device_add3s2(this%absvolflux_d, this%volflux1_d, &
974 this%volflux2_d, 1.0_rp, 1.0_rp, this%n_large)
975 call device_add2(this%absvolflux_d, this%volflux3_d, this%n_large)
976 call device_absval(this%absvolflux_d, this%n_large)
977 else
978 call pick_facet_value_hex(this%volflux1, u%x, &
979 this%lx, this%coef%msh%nelv)
980 call pick_facet_value_hex(this%volflux2, v%x, &
981 this%lx, this%coef%msh%nelv)
982 call pick_facet_value_hex(this%volflux3, w%x, &
983 this%lx, this%coef%msh%nelv)
984 call col2(this%volflux1, this%n1, this%n_large)
985 call col2(this%volflux2, this%n2, this%n_large)
986 call col2(this%volflux3, this%n3, this%n_large)
987 call add3(this%absvolflux, this%volflux1, this%volflux2, this%n_large)
988 call add2(this%absvolflux, this%volflux3, this%n_large)
989 call absval(this%absvolflux, this%n_large)
990 end if
991
992 end subroutine absvolflux_compute
993
999 subroutine pick_facet_value_hex(f_facet, f_field, lx, nelv)
1000 integer, intent(in) :: lx, nelv
1001 real(kind=rp), intent(in) :: f_field(lx, lx, lx, nelv)
1002 real(kind=rp), intent(inout) :: f_facet(lx + 2, lx + 2, lx + 2, nelv)
1003
1004 call copy(f_facet(1, 2: lx + 1, 2: lx + 1, :), &
1005 f_field(1, :, :, :), lx * lx * nelv)
1006 call copy(f_facet(lx + 2, 2: lx + 1, 2: lx + 1, :), &
1007 f_field(lx, :, :, :), lx * lx * nelv)
1008 call copy(f_facet(2: lx + 1, 1, 2: lx + 1, :), &
1009 f_field(:, 1, :, :), lx * lx * nelv)
1010 call copy(f_facet(2: lx + 1, lx + 2, 2: lx + 1, :), &
1011 f_field(:, lx, :, :), lx * lx * nelv)
1012 call copy(f_facet(2: lx + 1, 2: lx + 1, 1, :), &
1013 f_field(:, :, 1, :), lx * lx * nelv)
1014 call copy(f_facet(2: lx + 1, 2: lx + 1, lx + 2, :), &
1015 f_field(:, :, lx, :), lx * lx * nelv)
1016
1017 end subroutine pick_facet_value_hex
1018
1019end module gradient_jump_penalty
Map a Fortran array to a device (allocate and associate)
Definition device.F90:77
Copy data between host and device (or device and device)
Definition device.F90:71
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Coefficients.
Definition coef.f90:34
subroutine, public device_gradient_jump_penalty_finalize(penalty_d, penalty_facet_d, dphidxi_d, nx, nelv)
subroutine, public device_pick_facet_value_hex(b_d, a_d, nx, nelv)
subroutine, public device_add2s2(a_d, b_d, c1, n, strm)
Vector addition with scalar multiplication (multiplication on first argument)
subroutine, public device_add2(a_d, b_d, n, strm)
Vector addition .
subroutine, public device_add3s2(a_d, b_d, c_d, c1, c2, n, strm)
Returns .
subroutine, public device_col2(a_d, b_d, n, strm)
Vector multiplication .
subroutine, public device_absval(a_d, n, strm)
subroutine, public device_invcol2(a_d, b_d, n, strm)
Vector division .
subroutine, public device_col3(a_d, b_d, c_d, n, strm)
Vector multiplication with 3 vectors .
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:47
subroutine, public device_free(x_d)
Deallocate memory on the device.
Definition device.F90:219
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
Defines a field.
Definition field.f90:34
Gather-scatter.
Implements gradient_jump_penalty_t.
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...
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...
subroutine gradient_jump_penalty_compute_single(this, s)
Compute the gradient jump penalty term for a single field.
subroutine absvolflux_compute(this, u, v, w)
Compute the average of the volumetric flux over facets.
subroutine facet_factor_init(this)
Initialize the facet factor array.
subroutine g_compute(this, s)
Compute the average of the flux over facets.
real(kind=rp) function dist2_quadrature_hex(l, k, j, i, n, dm, coef)
subroutine pick_facet_value_hex(f_facet, f_field, lx, nelv)
Pick facet values of a field.
subroutine gradient_jump_penalty_init(this, json, fields, coef, variable_name)
Constructor.
subroutine eval_h2_hex(h2_el, n, i, coef)
Evaluate h^2 for each element for hexahedral mesh.
subroutine gradient_jump_penalty_compute(this, time)
Assign the gradient jump penalty term.
subroutine gradient_jump_penalty_init_from_components(this, coef, fields, start_time, end_time, a, b, variable_name)
Constructor from components.
subroutine gradient_jump_penalty_free(this)
Destructor for the gradient_jump_penalty_t class.
Defines Gather-scatter operations.
Definition gs_ops.f90:34
integer, parameter, public gs_op_add
Definition gs_ops.f90:36
Defines a hexahedron element.
Definition hex.f90:34
Utilities for retrieving parameters from the case files.
Definition math.f90:60
subroutine, public invcol2(a, b, n)
Vector division .
Definition math.f90:840
subroutine, public add3(a, b, c, n)
Vector addition .
Definition math.f90:739
subroutine, public add2(a, b, n)
Vector addition .
Definition math.f90:726
subroutine, public absval(a, n)
Take the absolute value of an array.
Definition math.f90:1372
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:854
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:249
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:867
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Operators.
Definition operators.f90:34
subroutine, public dudxyz(du, u, dr, ds, dt, coef)
Compute derivative of a scalar field along a single direction.
Definition operators.f90:90
Implements a point.
Definition point.f90:35
Defines a quadrilateral element.
Definition quad.f90:34
Defines a registry for storing solution fields.
Definition registry.f90:34
type(registry_t), target, public neko_registry
Global field registry.
Definition registry.f90:128
Implements the source_term_t type and a wrapper source_term_wrapper_t.
Defines a function space.
Definition space.f90:34
integer, parameter, public gll
Definition space.f90:49
Module with things related to the simulation time.
Utilities.
Definition utils.f90:35
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:56
Base type for an element.
Definition element.f90:44
field_list_t, To be able to group fields together
Hexahedron element.
Definition hex.f90:63
A point in with coordinates .
Definition point.f90:43
Quadrilateral element.
Definition quad.f90:58
Base abstract type for source terms.
The function space for the SEM solution fields.
Definition space.f90:63
A struct that contains all info about the time, expand as needed.