486 subroutine s_deltafunc(nBubs, lbk_rad, lbk_vel, lbk_s, updatedvar, kcomp)
488 integer,
intent(in) :: nBubs
489 real(wp),
dimension(1:lag_params%nBubs_glb,1:3,1:2),
intent(in) :: lbk_s
490 real(wp),
dimension(1:lag_params%nBubs_glb,1:2),
intent(in) :: lbk_rad, lbk_vel
491 type(scalar_field),
dimension(:),
intent(inout) :: updatedvar
492 type(scalar_field),
dimension(:),
intent(inout) :: kcomp
493 real(wp) :: strength_vel, strength_vol
494 real(wp) :: volpart, Vol
495 real(wp) :: y_kahan, t_kahan
496 integer :: i, j, k, lb, bub_idx
499# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
501# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
502#if defined(MFC_OpenACC)
503# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
505# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
506#elif defined(MFC_OpenMP)
507# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
509# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
511# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
513# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
515# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
517# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
519# 123 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
524 if (num_dims == 2)
then
525 vol = dx(i)*dy(j)*lag_params%charwidth
526 if (cyl_coord) vol = dx(i)*dy(j)*y_cc(j)*2._wp*pi
528 vol = dx(i)*dy(j)*dz(k)
533# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
534#if defined(MFC_OpenACC)
535# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
537# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
538#elif defined(MFC_OpenMP)
539# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
541# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
546 volpart = 4._wp/3._wp*pi*lbk_rad(bub_idx, 2)**3._wp
547 strength_vol = volpart
548 strength_vel = 4._wp*pi*lbk_rad(bub_idx, 2)**2._wp*lbk_vel(bub_idx, 2)
550 if (lag_params%kahan_summation)
then
552 y_kahan = real(strength_vol/vol, kind=wp) - kcomp(1)%sf(i, j, k)
553 t_kahan = updatedvar(1)%sf(i, j, k) + y_kahan
554 kcomp(1)%sf(i, j, k) = (t_kahan - updatedvar(1)%sf(i, j, k)) - y_kahan
555 updatedvar(1)%sf(i, j, k) = t_kahan
558 y_kahan = real(strength_vel/vol, kind=wp) - kcomp(2)%sf(i, j, k)
559 t_kahan = updatedvar(2)%sf(i, j, k) + y_kahan
560 kcomp(2)%sf(i, j, k) = (t_kahan - updatedvar(2)%sf(i, j, k)) - y_kahan
561 updatedvar(2)%sf(i, j, k) = t_kahan
563 updatedvar(1)%sf(i, j, k) = updatedvar(1)%sf(i, j, k) + real(strength_vol/vol, kind=wp)
564 updatedvar(2)%sf(i, j, k) = updatedvar(2)%sf(i, j, k) + real(strength_vel/vol, kind=wp)
568 if (lag_params%kahan_summation .and. lag_params%cluster_type >= 4)
then
569 y_kahan = real((strength_vol*strength_vel)/vol, kind=wp) - kcomp(5)%sf(i, j, k)
570 t_kahan = updatedvar(5)%sf(i, j, k) + y_kahan
571 kcomp(5)%sf(i, j, k) = (t_kahan - updatedvar(5)%sf(i, j, k)) - y_kahan
572 updatedvar(5)%sf(i, j, k) = t_kahan
573 else if (lag_params%cluster_type >= 4)
then
574 updatedvar(5)%sf(i, j, k) = updatedvar(5)%sf(i, j, k) + real((strength_vol*strength_vel)/vol, kind=wp)
581# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
582#if defined(MFC_OpenACC)
583# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
585# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
586#elif defined(MFC_OpenMP)
587# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
589# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
591# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
598 subroutine s_gaussian(nBubs, lbk_rad, lbk_vel, lbk_s, lbk_pos, updatedvar, kcomp)
600 integer,
intent(in) :: nBubs
601 real(wp),
dimension(1:lag_params%nBubs_glb,1:3,1:2),
intent(in) :: lbk_s, lbk_pos
602 real(wp),
dimension(1:lag_params%nBubs_glb,1:2),
intent(in) :: lbk_rad, lbk_vel
603 type(scalar_field),
dimension(:),
intent(inout) :: updatedvar
604 type(scalar_field),
dimension(:),
intent(inout) :: kcomp
605 real(wp),
dimension(3) :: center, nodecoord, s_coord
606 integer,
dimension(3) :: cell, cellijk
607 real(wp) :: stddsv, volpart
608 real(wp) :: strength_vel, strength_vol
609 real(wp) :: func, func2
610 real(wp) :: y_kahan, t_kahan
611 integer :: i, j, k, di, dj, dk, lb, bub_idx
612 integer :: di_beg, di_end, dj_beg, dj_end, dk_beg, dk_end
613 integer :: smear_x_beg, smear_x_end
614 integer :: smear_y_beg, smear_y_end
615 integer :: smear_z_beg, smear_z_end
619 smear_x_beg = -mapcells - 1
620 smear_x_end = m + mapcells + 1
621 smear_y_beg = merge(-mapcells - 1, 0, n > 0)
622 smear_y_end = merge(n + mapcells + 1, n, n > 0)
623 smear_z_beg = merge(-mapcells - 1, 0, p > 0)
624 smear_z_end = merge(p + mapcells + 1, p, p > 0)
627# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
629# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
630#if defined(MFC_OpenACC)
631# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
633# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
635# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
636#elif defined(MFC_OpenMP)
637# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
639# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
641# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
643# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
645# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
647# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
649# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
651# 210 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
652 do k = smear_z_beg, smear_z_end
653 do j = smear_y_beg, smear_y_end
654 do i = smear_x_beg, smear_x_end
659 nodecoord(1) = x_cc(i)
660 nodecoord(2) = y_cc(j)
662 if (p > 0) nodecoord(3) = z_cc(k)
665 di_beg = max(i - mapcells, 0)
666 di_end = min(i + mapcells, m)
667 dj_beg = max(j - mapcells, 0)
668 dj_end = min(j + mapcells, n)
669 dk_beg = max(k - mapcells, 0)
670 dk_end = min(k + mapcells, p)
673# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
674#if defined(MFC_OpenACC)
675# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
677# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
678#elif defined(MFC_OpenMP)
679# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
681# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
683 do dk = dk_beg, dk_end
685# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
686#if defined(MFC_OpenACC)
687# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
689# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
690#elif defined(MFC_OpenMP)
691# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
693# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
695 do dj = dj_beg, dj_end
697# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
698#if defined(MFC_OpenACC)
699# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
701# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
702#elif defined(MFC_OpenMP)
703# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
705# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
707 do di = di_beg, di_end
709# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
710#if defined(MFC_OpenACC)
711# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
713# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
714#elif defined(MFC_OpenMP)
715# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
717# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
723 volpart = 4._wp/3._wp*pi*lbk_rad(bub_idx, 2)**3._wp
724 s_coord(1:3) = lbk_s(bub_idx,1:3,2)
728 strength_vol = volpart
729 strength_vel = 4._wp*pi*lbk_rad(bub_idx, 2)**2._wp*lbk_vel(bub_idx, 2)
731 center(1:2) = lbk_pos(bub_idx,1:2,2)
733 if (p > 0) center(3) = lbk_pos(bub_idx, 3, 2)
738 y_kahan = real(func*strength_vol, kind=wp) - kcomp(1)%sf(i, j, k)
739 t_kahan = updatedvar(1)%sf(i, j, k) + y_kahan
740 kcomp(1)%sf(i, j, k) = (t_kahan - updatedvar(1)%sf(i, j, k)) - y_kahan
741 updatedvar(1)%sf(i, j, k) = t_kahan
744 y_kahan = real(func*strength_vel, kind=wp) - kcomp(2)%sf(i, j, k)
745 t_kahan = updatedvar(2)%sf(i, j, k) + y_kahan
746 kcomp(2)%sf(i, j, k) = (t_kahan - updatedvar(2)%sf(i, j, k)) - y_kahan
747 updatedvar(2)%sf(i, j, k) = t_kahan
749 if (lag_params%cluster_type >= 4)
then
751 y_kahan = real(func2*strength_vol*strength_vel, kind=wp) - kcomp(5)%sf(i, j, k)
752 t_kahan = updatedvar(5)%sf(i, j, k) + y_kahan
753 kcomp(5)%sf(i, j, k) = (t_kahan - updatedvar(5)%sf(i, j, k)) - y_kahan
754 updatedvar(5)%sf(i, j, k) = t_kahan
764# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
765#if defined(MFC_OpenACC)
766# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
768# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
769#elif defined(MFC_OpenMP)
770# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
772# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
774# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
783# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
785# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
787# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
789# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
791# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
793# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
795# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
797# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
799# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
801# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
803# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
805# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
807# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
809# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
811# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
813# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
815# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
818 real(wp),
dimension(3),
intent(in) :: center
819 integer,
dimension(3),
intent(in) :: cellaux
820 real(wp),
dimension(3),
intent(in) :: nodecoord
821 real(wp),
intent(in) :: stddsv
822 real(wp),
intent(in) :: strength_idx
823 real(wp),
intent(out) :: func
826 real(wp) :: theta, dtheta, L2, dzp, Lz2, zc
827 real(wp) :: Nr, Nr_count
829 distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + (center(3) - nodecoord(3))**2._wp)
831 if (num_dims == 3)
then
833 func = exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
839 nr = ceiling(2._wp*pi*nodecoord(2)/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
841 l2 = center(2)**2._wp + nodecoord(2)**2._wp - 2._wp*center(2)*nodecoord(2)*cos(theta)
842 distance = sqrt((center(1) - nodecoord(1))**2._wp + l2)
844 func = dtheta/2._wp/pi*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
846 do while (nr_count < nr - 1._wp)
847 nr_count = nr_count + 1._wp
848 theta = nr_count*dtheta
850 l2 = center(2)**2._wp + nodecoord(2)**2._wp - 2._wp*center(2)*nodecoord(2)*cos(theta)
851 distance = sqrt((center(1) - nodecoord(1))**2._wp + l2)
853 func = func + dtheta/2._wp/pi*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv) &
854 & **(3._wp*(strength_idx + 1._wp))
859 dzp = (lag_params%charwidth/(lag_params%charNz + 1._wp))
862 do i = 0, lag_params%charNz
863 zc = (-lag_params%charwidth/2._wp + dzp*(0.5_wp + i))
864 lz2 = (center(3) - zc)**2._wp
865 distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + lz2)
866 func = func + dzp/lag_params%charwidth*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
1365# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1367# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1369# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1371# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1373# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1375# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1377# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1379# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1381# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1383# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1385# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1387# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1389# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1391# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1393# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1395# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1397# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1400 real(wp),
intent(in) :: pos
1401 integer,
dimension(3),
intent(in) :: cell
1402 integer,
intent(in) :: i
1403 type(scalar_field),
dimension(sys_size),
intent(in) :: q_prim_vf
1405 real(wp),
dimension(5) :: xi, eta,
l
1407 if (fd_order == 2)
then
1409 xi(1) = x_cc(cell(1) - 1)
1410 eta(1) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 1, cell(2), cell(3))
1411 xi(2) = x_cc(cell(1))
1412 eta(2) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1), cell(2), cell(3))
1413 xi(3) = x_cc(cell(1) + 1)
1414 eta(3) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 1, cell(2), cell(3))
1415 else if (i == 2)
then
1416 xi(1) = y_cc(cell(2) - 1)
1417 eta(1) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 1, cell(3))
1418 xi(2) = y_cc(cell(2))
1419 eta(2) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2), cell(3))
1420 xi(3) = y_cc(cell(2) + 1)
1421 eta(3) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 1, cell(3))
1422 else if (i == 3)
then
1423 xi(1) = z_cc(cell(3) - 1)
1424 eta(1) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 1)
1425 xi(2) = z_cc(cell(3))
1426 eta(2) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3))
1427 xi(3) = z_cc(cell(3) + 1)
1428 eta(3) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 1)
1431 l(1) = ((pos - xi(2))*(pos - xi(3)))/((xi(1) - xi(2))*(xi(1) - xi(3)))
1432 l(2) = ((pos - xi(1))*(pos - xi(3)))/((xi(2) - xi(1))*(xi(2) - xi(3)))
1433 l(3) = ((pos - xi(1))*(pos - xi(2)))/((xi(3) - xi(1))*(xi(3) - xi(2)))
1435 v =
l(1)*eta(1) +
l(2)*eta(2) +
l(3)*eta(3)
1436 else if (fd_order == 4)
then
1438 xi(1) = x_cc(cell(1) - 2)
1439 eta(1) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 2, cell(2), cell(3))
1440 xi(2) = x_cc(cell(1) - 1)
1441 eta(2) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 1, cell(2), cell(3))
1442 xi(3) = x_cc(cell(1))
1443 eta(3) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1), cell(2), cell(3))
1444 xi(4) = x_cc(cell(1) + 1)
1445 eta(4) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 1, cell(2), cell(3))
1446 xi(5) = x_cc(cell(1) + 2)
1447 eta(5) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 2, cell(2), cell(3))
1448 else if (i == 2)
then
1449 xi(1) = y_cc(cell(2) - 2)
1450 eta(1) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 2, cell(3))
1451 xi(2) = y_cc(cell(2) - 1)
1452 eta(2) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 1, cell(3))
1453 xi(3) = y_cc(cell(2))
1454 eta(3) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2), cell(3))
1455 xi(4) = y_cc(cell(2) + 1)
1456 eta(4) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 1, cell(3))
1457 xi(5) = y_cc(cell(2) + 2)
1458 eta(5) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 2, cell(3))
1459 else if (i == 3)
then
1460 xi(1) = z_cc(cell(3) - 2)
1461 eta(1) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 2)
1462 xi(2) = z_cc(cell(3) - 1)
1463 eta(2) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 1)
1464 xi(3) = z_cc(cell(3))
1465 eta(3) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3))
1466 xi(4) = z_cc(cell(3) + 1)
1467 eta(4) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 1)
1468 xi(5) = z_cc(cell(3) + 2)
1469 eta(5) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 2)
1472 l(1) = ((pos - xi(2))*(pos - xi(3))*(pos - xi(4))*(pos - xi(5)))/((xi(1) - xi(2))*(xi(1) - xi(3))*(xi(1) - xi(4)) &
1474 l(2) = ((pos - xi(1))*(pos - xi(3))*(pos - xi(4))*(pos - xi(5)))/((xi(2) - xi(1))*(xi(2) - xi(3))*(xi(2) - xi(4)) &
1476 l(3) = ((pos - xi(1))*(pos - xi(2))*(pos - xi(4))*(pos - xi(5)))/((xi(3) - xi(1))*(xi(3) - xi(2))*(xi(3) - xi(4)) &
1478 l(4) = ((pos - xi(1))*(pos - xi(2))*(pos - xi(3))*(pos - xi(5)))/((xi(4) - xi(1))*(xi(4) - xi(2))*(xi(4) - xi(3)) &
1480 l(5) = ((pos - xi(1))*(pos - xi(2))*(pos - xi(3))*(pos - xi(4)))/((xi(5) - xi(1))*(xi(5) - xi(2))*(xi(5) - xi(3)) &
1483 v =
l(1)*eta(1) +
l(2)*eta(2) +
l(3)*eta(3) +
l(4)*eta(4) +
l(5)*eta(5)
1493 function f_get_bubble_force(pos, rad, rdot, vel, mg, mv, Re, rho, cell, i, q_prim_vf)
result(force)
1496# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1498# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1500# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1502# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1504# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1506# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1508# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1510 real(wp),
intent(in) :: pos, rad, rdot, mg, mv, re, rho, vel
1511 integer,
dimension(3),
intent(in) :: cell
1512 integer,
intent(in) :: i
1513 type(scalar_field),
dimension(sys_size),
intent(in) :: q_prim_vf
1514 real(wp) :: dp, vol, force
1517 if (fd_order > 1)
then
1520 v_rel = vel - q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(cell(1), cell(2), cell(3))
1525 if (lag_params%drag_model == 1)
then
1526 force = force - (4._wp*pi*rad*v_rel)/re
1527 else if (lag_params%drag_model == 2)
then
1528 force = force - (6._wp*pi*rad*v_rel)/re
1529 else if (lag_params%drag_model == 3)
then
1530 force = force - (12._wp*pi*rad*v_rel)/re
1533 if (lag_pressure_force)
then
1536 dp =
grad_p_x(cell(1), cell(2), cell(3))
1537 else if (i == 2)
then
1538 dp =
grad_p_y(cell(1), cell(2), cell(3))
1539 else if (i == 3)
then
1540 dp =
grad_p_z(cell(1), cell(2), cell(3))
1543 vol = (4._wp/3._wp)*pi*(rad**3._wp)
1544 force = force - vol*dp
1547 if (lag_params%gravity_force)
then
1548 force = force + (mg + mv)*accel_bf(i)