Neko 1.99.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!
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
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_field_registry%get_field("u")
215 this%v => neko_field_registry%get_field("v")
216 this%w => neko_field_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_field_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
795 end subroutine gradient_jump_penalty_free
796
800 class(gradient_jump_penalty_t), intent(inout) :: this
801 type(field_t), intent(in) :: s
802
803 class(element_t), pointer :: ep
804 integer :: i
805
806 call g_compute(this, s)
807 call absvolflux_compute(this, this%u, this%v, this%w)
808
809 if (neko_bcknd_device .eq. 1) then
810 call device_col3(this%penalty_facet_d, this%absvolflux_d, this%G_d, &
811 this%n_large)
812 call device_col2(this%penalty_facet_d, this%facet_factor_d, this%n_large)
813 call device_gradient_jump_penalty_finalize(this%penalty_d, &
814 this%penalty_facet_d, &
815 this%dphidxi_d, &
816 this%lx, this%coef%msh%nelv)
817 else
818 call col3(this%penalty_facet, this%absvolflux, this%G, this%n_large)
819 call col2(this%penalty_facet, this%facet_factor, this%n_large)
820 call gradient_jump_penalty_finalize(this%penalty, this%penalty_facet, &
821 this%dphidxi, &
822 this%lx, this%coef%msh%nelv)
823 end if
824
826
829 subroutine gradient_jump_penalty_compute(this, time)
830 class(gradient_jump_penalty_t), intent(inout) :: this
831 type(time_state_t), intent(in) :: time
832 integer :: i, n_fields, n
833
834 n_fields = this%fields%size()
835 n = this%coef%dof%size()
836
837
838
839 do i = 1, n_fields
840
841 call this%compute_single(this%s_fields%items(i)%ptr)
842
843 if (neko_bcknd_device .eq. 1) then
844 call device_invcol2(this%penalty_d, this%coef%B_d, n)
845 call device_add2(this%fields%x_d(i), this%penalty_d, n)
846 else
847 call invcol2(this%penalty, this%coef%B, n)
848 call add2(this%fields%items(i)%ptr%x, this%penalty, n)
849 end if
850 end do
851
852 end subroutine gradient_jump_penalty_compute
853
861 subroutine gradient_jump_penalty_finalize(penalty, wa, dphidxi, lx, nelv)
862 integer, intent(in) :: lx, nelv
863 real(kind=rp), intent(inout) :: penalty(lx, lx, lx, nelv)
864 real(kind=rp), intent(in) :: wa(lx + 2, lx + 2, lx + 2, nelv)
865 real(kind=rp), intent(in) :: dphidxi(lx, lx)
866
867 call gradient_jump_penalty_finalize_hex(penalty, wa, dphidxi, lx, nelv)
868
869 end subroutine gradient_jump_penalty_finalize
870
878 subroutine gradient_jump_penalty_finalize_hex(penalty, wa, dphidxi, lx, nelv)
879 integer, intent(in) :: lx, nelv
880 real(kind=rp), intent(inout) :: penalty(lx, lx, lx, nelv)
881 real(kind=rp), intent(in) :: wa(lx + 2, lx + 2, lx + 2, nelv)
882 real(kind=rp), intent(in) :: dphidxi(lx, lx)
883
884 integer :: i, j, k
885
886 do i = 1, lx
887 do j = 1, lx
888 do k = 1, lx
889 penalty(i, j, k, :) = &
890 wa(1, j + 1, k + 1, :) * &
891 dphidxi(1, i) + &
892 wa(lx + 2, j + 1, k + 1, :) * &
893 dphidxi(lx, i) + &
894 wa(i + 1, 1, k + 1, :) * &
895 dphidxi(1, j) + &
896 wa(i + 1, lx + 2, k + 1, :) * &
897 dphidxi(lx, j) + &
898 wa(i + 1, j + 1, 1, :) * &
899 dphidxi(1, k) + &
900 wa(i + 1, j + 1, lx + 2, :) * &
901 dphidxi(lx, k)
902 end do
903 end do
904 end do
905
907
910 subroutine g_compute(this, s)
911 class(gradient_jump_penalty_t), intent(inout) :: this
912 type(field_t), intent(in) :: s
913
914 call dudxyz(this%grad1, s%x, this%coef%drdx, &
915 this%coef%dsdx, this%coef%dtdx, this%coef)
916 call dudxyz(this%grad2, s%x, this%coef%drdy, &
917 this%coef%dsdy, this%coef%dtdy, this%coef)
918 call dudxyz(this%grad3, s%x, this%coef%drdz, &
919 this%coef%dsdz, this%coef%dtdz, this%coef)
920
921 if (neko_bcknd_device .eq. 1) then
922 call device_pick_facet_value_hex(this%flux1_d, this%grad1_d, &
923 this%lx, this%coef%msh%nelv)
924 call device_pick_facet_value_hex(this%flux2_d, this%grad2_d, &
925 this%lx, this%coef%msh%nelv)
926 call device_pick_facet_value_hex(this%flux3_d, this%grad3_d, &
927 this%lx, this%coef%msh%nelv)
928 call device_col2(this%flux1_d, this%n1_d, this%n_large)
929 call device_col2(this%flux2_d, this%n2_d, this%n_large)
930 call device_col2(this%flux3_d, this%n3_d, this%n_large)
931 call device_add3s2(this%G_d, this%flux1_d, this%flux2_d, &
932 1.0_rp, 1.0_rp, this%n_large)
933 call device_add2(this%G_d, this%flux3_d, this%n_large)
934 else
935 call pick_facet_value_hex(this%flux1, this%grad1, &
936 this%lx, this%coef%msh%nelv)
937 call pick_facet_value_hex(this%flux2, this%grad2, &
938 this%lx, this%coef%msh%nelv)
939 call pick_facet_value_hex(this%flux3, this%grad3, &
940 this%lx, this%coef%msh%nelv)
941 call col2(this%flux1, this%n1, this%n_large)
942 call col2(this%flux2, this%n2, this%n_large)
943 call col2(this%flux3, this%n3, this%n_large)
944 call add3(this%G, this%flux1, this%flux2, this%n_large)
945 call add2(this%G, this%flux3, this%n_large)
946 end if
947
948 call this%gs_GJP%op(this%G, this%n_large, gs_op_add)
949
950 end subroutine g_compute
951
956 subroutine absvolflux_compute(this, u, v, w)
957 class(gradient_jump_penalty_t), intent(inout) :: this
958 type(field_t), intent(in) :: u, v, w
959
960 integer :: i
961
962 if (neko_bcknd_device .eq. 1) then
963 call device_pick_facet_value_hex(this%volflux1_d, u%x_d, this%lx, &
964 this%coef%msh%nelv)
965 call device_pick_facet_value_hex(this%volflux2_d, v%x_d, this%lx, &
966 this%coef%msh%nelv)
967 call device_pick_facet_value_hex(this%volflux3_d, w%x_d, this%lx, &
968 this%coef%msh%nelv)
969 call device_col2(this%volflux1_d, this%n1_d, this%n_large)
970 call device_col2(this%volflux2_d, this%n2_d, this%n_large)
971 call device_col2(this%volflux3_d, this%n3_d, this%n_large)
972 call device_add3s2(this%absvolflux_d, this%volflux1_d, &
973 this%volflux2_d, 1.0_rp, 1.0_rp, this%n_large)
974 call device_add2(this%absvolflux_d, this%volflux3_d, this%n_large)
975 call device_absval(this%absvolflux_d, this%n_large)
976 else
977 call pick_facet_value_hex(this%volflux1, u%x, &
978 this%lx, this%coef%msh%nelv)
979 call pick_facet_value_hex(this%volflux2, v%x, &
980 this%lx, this%coef%msh%nelv)
981 call pick_facet_value_hex(this%volflux3, w%x, &
982 this%lx, this%coef%msh%nelv)
983 call col2(this%volflux1, this%n1, this%n_large)
984 call col2(this%volflux2, this%n2, this%n_large)
985 call col2(this%volflux3, this%n3, this%n_large)
986 call add3(this%absvolflux, this%volflux1, this%volflux2, this%n_large)
987 call add2(this%absvolflux, this%volflux3, this%n_large)
988 call absval(this%absvolflux, this%n_large)
989 end if
990
991 end subroutine absvolflux_compute
992
998 subroutine pick_facet_value_hex(f_facet, f_field, lx, nelv)
999 integer, intent(in) :: lx, nelv
1000 real(kind=rp), intent(in) :: f_field(lx, lx, lx, nelv)
1001 real(kind=rp), intent(inout) :: f_facet(lx + 2, lx + 2, lx + 2, nelv)
1002
1003 call copy(f_facet(1, 2: lx + 1, 2: lx + 1, :), &
1004 f_field(1, :, :, :), lx * lx * nelv)
1005 call copy(f_facet(lx + 2, 2: lx + 1, 2: lx + 1, :), &
1006 f_field(lx, :, :, :), lx * lx * nelv)
1007 call copy(f_facet(2: lx + 1, 1, 2: lx + 1, :), &
1008 f_field(:, 1, :, :), lx * lx * nelv)
1009 call copy(f_facet(2: lx + 1, lx + 2, 2: lx + 1, :), &
1010 f_field(:, lx, :, :), lx * lx * nelv)
1011 call copy(f_facet(2: lx + 1, 2: lx + 1, 1, :), &
1012 f_field(:, :, 1, :), lx * lx * nelv)
1013 call copy(f_facet(2: lx + 1, 2: lx + 1, lx + 2, :), &
1014 f_field(:, :, lx, :), lx * lx * nelv)
1015
1016 end subroutine pick_facet_value_hex
1017
1018end module gradient_jump_penalty
Map a Fortran array to a device (allocate and associate)
Definition device.F90:72
Copy data between host and device (or device and device)
Definition device.F90:66
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:214
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
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:846
subroutine, public add3(a, b, c, n)
Vector addition .
Definition math.f90:745
subroutine, public add2(a, b, n)
Vector addition .
Definition math.f90:732
subroutine, public absval(a, n)
Take the absolute value of an array.
Definition math.f90:1331
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:860
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:255
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:873
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:80
Implements a point.
Definition point.f90:35
Defines a quadrilateral element.
Definition quad.f90:34
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:48
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:55
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:62
A struct that contains all info about the time, expand as needed.