197 start_time, end_time, a, b, variable_name)
199 type(
coef_t),
target,
intent(in) :: coef
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
206 integer :: i, j, k, l
207 real(kind=
rp),
allocatable :: zg(:)
208 real(kind=
rp) :: normal(3)
212 call this%init_base(fields, coef, start_time, end_time)
218 if (fields%size() .eq. 1)
then
219 call this%s_fields%init(1)
220 call this%s_fields%assign(1, &
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)
228 call neko_error(
"The GJP source assumes either 3 or 1 RHS fields.")
234 this%p = coef%dof%Xh%lx - 1
235 this%lx = coef%dof%Xh%lx
237 if (this%p .gt. 1)
then
238 this%tau = -a * (this%p + 1) ** (-b)
243 this%n = this%lx ** 3 * this%coef%msh%nelv
244 this%n_large = (this%lx + 2) ** 3 * this%coef%msh%nelv
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)
253 &supported now for gradient jump penalty")
256 this%n_facet_max = maxval(this%n_facet)
258 allocate(this%h2(this%lx + 2, this%lx + 2, &
259 this%lx + 2, this%coef%msh%nelv))
261 do i = 1, this%coef%msh%nelv
262 select type (ep => this%coef%msh%elements(i)%e)
264 call eval_h2_hex(this%h2(:, :, :, i), this%lx, i, this%coef)
266 call neko_error(
"Gradient jump penalty error: mesh size &
267 &evaluation is not supported for quad_t")
271 allocate(zg(this%lx))
272 allocate(this%dphidxi(this%lx, this%lx))
277 this%dphidxi(j,i) = this%coef%Xh%dx(j,i)
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))
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))
312 do i = 1, this%coef%msh%nelv
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)
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)
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)
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)
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)
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)
348 call neko_error(
"The face index is not correct")
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)
365 call device_map(this%dphidxi, this%dphidxi_d, &
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)
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)
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)
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)
397 call device_memcpy(this%facet_factor, this%facet_factor_d, this%n_large,&
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
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
464 normal1 = coef%get_normal(1, l, k, i, 1)
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)
475 x2 = dm%x(n, l, k, i)
476 y2 = dm%y(n, l, k, i)
477 z2 = dm%z(n, l, k, i)
481 normal1 = coef%get_normal(1, l, k, i, 2)
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)
492 x2 = dm%x(1, l, k, i)
493 y2 = dm%y(1, l, k, i)
494 z2 = dm%z(1, l, k, i)
496 normal1 = coef%get_normal(1, l, k, i, 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)
507 x2 = dm%x(l, n, k, i)
508 y2 = dm%y(l, n, k, i)
509 z2 = dm%z(l, n, k, i)
511 normal1 = coef%get_normal(1, l, k, i, 4)
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)
522 x2 = dm%x(l, 1, k, i)
523 y2 = dm%y(l, 1, k, i)
524 z2 = dm%z(l, 1, k, i)
526 normal1 = coef%get_normal(1, l, k, i, 5)
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)
537 x2 = dm%x(l, k, n, i)
538 y2 = dm%y(l, k, n, i)
539 z2 = dm%z(l, k, n, i)
541 normal1 = coef%get_normal(1, l, k, i, 6)
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)
552 x2 = dm%x(l, k, 1, i)
553 y2 = dm%y(l, k, 1, i)
554 z2 = dm%z(l, k, 1, i)
556 call neko_error(
"The face index is not correct")
564 dist_1 = v1*n11 + v2*n12 + v3*n13
565 dist_2 = - (v1*n21 + v2*n22 + v3*n23)
567 dist2 = ((dist_1 + dist_2)/2.0_rp)*((dist_1 + dist_2)/2.0_rp)
574 integer :: i, j, k, l
575 real(kind=
rp) :: area_tmp
577 allocate(this%facet_factor(this%lx + 2, this%lx + 2, &
578 this%lx + 2, this%coef%msh%nelv))
580 associate(facet_factor => this%facet_factor, &
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)
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) ) &
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)
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) ) &
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)
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) ) &
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) * &
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) &
651 * jacinv(l, k, lx, i)
653 call neko_error(
"The face index is not correct")
912 type(field_t),
intent(in) :: s
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)
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)
936 this%lx, this%coef%msh%nelv)
938 this%lx, this%coef%msh%nelv)
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)
948 call this%gs_GJP%op(this%G, this%n_large, gs_op_add)
958 type(field_t),
intent(in) :: u, v, w
962 if (neko_bcknd_device .eq. 1)
then
963 call device_pick_facet_value_hex(this%volflux1_d, u%x_d, this%lx, &
965 call device_pick_facet_value_hex(this%volflux2_d, v%x_d, this%lx, &
967 call device_pick_facet_value_hex(this%volflux3_d, w%x_d, this%lx, &
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)
978 this%lx, this%coef%msh%nelv)
980 this%lx, this%coef%msh%nelv)
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)
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)
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)