474 subroutine s_deltafunc(nBubs, lbk_rad, lbk_vel, lbk_s, updatedvar, kcomp)
476 integer,
intent(in) :: nBubs
477 real(wp),
dimension(1:lag_params%nBubs_glb,1:3,1:2),
intent(in) :: lbk_s
478 real(wp),
dimension(1:lag_params%nBubs_glb,1:2),
intent(in) :: lbk_rad, lbk_vel
479 type(scalar_field),
dimension(:),
intent(inout) :: updatedvar
480 type(scalar_field),
dimension(:),
intent(inout) :: kcomp
481 real(wp) :: strength_vel, strength_vol
482 real(wp) :: volpart, Vol
483 real(wp) :: y_kahan, t_kahan
484 integer :: i, j, k, lb, bub_idx
487# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
489# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
490#if defined(MFC_OpenACC)
491# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
493# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
494#elif defined(MFC_OpenMP)
495# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
497# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
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"
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"
507# 123 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
512 if (num_dims == 2)
then
513 vol = dx(i)*dy(j)*lag_params%charwidth
514 if (cyl_coord) vol = dx(i)*dy(j)*y_cc(j)*2._wp*pi
516 vol = dx(i)*dy(j)*dz(k)
521# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
522#if defined(MFC_OpenACC)
523# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
525# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
526#elif defined(MFC_OpenMP)
527# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
529# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
534 volpart = 4._wp/3._wp*pi*lbk_rad(bub_idx, 2)**3._wp
535 strength_vol = volpart
536 strength_vel = 4._wp*pi*lbk_rad(bub_idx, 2)**2._wp*lbk_vel(bub_idx, 2)
538 if (lag_params%kahan_summation)
then
540 y_kahan = real(strength_vol/vol, kind=wp) - kcomp(1)%sf(i, j, k)
541 t_kahan = updatedvar(1)%sf(i, j, k) + y_kahan
542 kcomp(1)%sf(i, j, k) = (t_kahan - updatedvar(1)%sf(i, j, k)) - y_kahan
543 updatedvar(1)%sf(i, j, k) = t_kahan
546 y_kahan = real(strength_vel/vol, kind=wp) - kcomp(2)%sf(i, j, k)
547 t_kahan = updatedvar(2)%sf(i, j, k) + y_kahan
548 kcomp(2)%sf(i, j, k) = (t_kahan - updatedvar(2)%sf(i, j, k)) - y_kahan
549 updatedvar(2)%sf(i, j, k) = t_kahan
551 updatedvar(1)%sf(i, j, k) = updatedvar(1)%sf(i, j, k) + real(strength_vol/vol, kind=wp)
552 updatedvar(2)%sf(i, j, k) = updatedvar(2)%sf(i, j, k) + real(strength_vel/vol, kind=wp)
556 if (lag_params%kahan_summation .and. lag_params%cluster_type >= 4)
then
557 y_kahan = real((strength_vol*strength_vel)/vol, kind=wp) - kcomp(5)%sf(i, j, k)
558 t_kahan = updatedvar(5)%sf(i, j, k) + y_kahan
559 kcomp(5)%sf(i, j, k) = (t_kahan - updatedvar(5)%sf(i, j, k)) - y_kahan
560 updatedvar(5)%sf(i, j, k) = t_kahan
561 else if (lag_params%cluster_type >= 4)
then
562 updatedvar(5)%sf(i, j, k) = updatedvar(5)%sf(i, j, k) + real((strength_vol*strength_vel)/vol, kind=wp)
569# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
570#if defined(MFC_OpenACC)
571# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
573# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
574#elif defined(MFC_OpenMP)
575# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
577# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
579# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
586 subroutine s_gaussian(nBubs, lbk_rad, lbk_vel, lbk_s, lbk_pos, updatedvar, kcomp)
588 integer,
intent(in) :: nBubs
589 real(wp),
dimension(1:lag_params%nBubs_glb,1:3,1:2),
intent(in) :: lbk_s, lbk_pos
590 real(wp),
dimension(1:lag_params%nBubs_glb,1:2),
intent(in) :: lbk_rad, lbk_vel
591 type(scalar_field),
dimension(:),
intent(inout) :: updatedvar
592 type(scalar_field),
dimension(:),
intent(inout) :: kcomp
593 real(wp),
dimension(3) :: center, nodecoord, s_coord
594 integer,
dimension(3) :: cell, cellijk
595 real(wp) :: stddsv, volpart
596 real(wp) :: strength_vel, strength_vol
597 real(wp) :: func, func2
598 real(wp) :: y_kahan, t_kahan
599 integer :: i, j, k, di, dj, dk, lb, bub_idx
600 integer :: di_beg, di_end, dj_beg, dj_end, dk_beg, dk_end
601 integer :: smear_x_beg, smear_x_end
602 integer :: smear_y_beg, smear_y_end
603 integer :: smear_z_beg, smear_z_end
607 smear_x_beg = -mapcells - 1
608 smear_x_end = m + mapcells + 1
609 smear_y_beg = merge(-mapcells - 1, 0, n > 0)
610 smear_y_end = merge(n + mapcells + 1, n, n > 0)
611 smear_z_beg = merge(-mapcells - 1, 0, p > 0)
612 smear_z_end = merge(p + mapcells + 1, p, p > 0)
615# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
617# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
618#if defined(MFC_OpenACC)
619# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
621# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
623# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
624#elif defined(MFC_OpenMP)
625# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
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"
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"
637# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
639# 210 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
640 do k = smear_z_beg, smear_z_end
641 do j = smear_y_beg, smear_y_end
642 do i = smear_x_beg, smear_x_end
647 nodecoord(1) = x_cc(i)
648 nodecoord(2) = y_cc(j)
650 if (p > 0) nodecoord(3) = z_cc(k)
653 di_beg = max(i - mapcells, 0)
654 di_end = min(i + mapcells, m)
655 dj_beg = max(j - mapcells, 0)
656 dj_end = min(j + mapcells, n)
657 dk_beg = max(k - mapcells, 0)
658 dk_end = min(k + mapcells, p)
661# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
662#if defined(MFC_OpenACC)
663# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
665# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
666#elif defined(MFC_OpenMP)
667# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
669# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
671 do dk = dk_beg, dk_end
673# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
674#if defined(MFC_OpenACC)
675# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
677# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
678#elif defined(MFC_OpenMP)
679# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
681# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
683 do dj = dj_beg, dj_end
685# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
686#if defined(MFC_OpenACC)
687# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
689# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
690#elif defined(MFC_OpenMP)
691# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
693# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
695 do di = di_beg, di_end
697# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
698#if defined(MFC_OpenACC)
699# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
701# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
702#elif defined(MFC_OpenMP)
703# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
705# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
711 volpart = 4._wp/3._wp*pi*lbk_rad(bub_idx, 2)**3._wp
712 s_coord(1:3) = lbk_s(bub_idx,1:3,2)
716 strength_vol = volpart
717 strength_vel = 4._wp*pi*lbk_rad(bub_idx, 2)**2._wp*lbk_vel(bub_idx, 2)
719 center(1:2) = lbk_pos(bub_idx,1:2,2)
721 if (p > 0) center(3) = lbk_pos(bub_idx, 3, 2)
726 y_kahan = real(func*strength_vol, kind=wp) - kcomp(1)%sf(i, j, k)
727 t_kahan = updatedvar(1)%sf(i, j, k) + y_kahan
728 kcomp(1)%sf(i, j, k) = (t_kahan - updatedvar(1)%sf(i, j, k)) - y_kahan
729 updatedvar(1)%sf(i, j, k) = t_kahan
732 y_kahan = real(func*strength_vel, kind=wp) - kcomp(2)%sf(i, j, k)
733 t_kahan = updatedvar(2)%sf(i, j, k) + y_kahan
734 kcomp(2)%sf(i, j, k) = (t_kahan - updatedvar(2)%sf(i, j, k)) - y_kahan
735 updatedvar(2)%sf(i, j, k) = t_kahan
737 if (lag_params%cluster_type >= 4)
then
739 y_kahan = real(func2*strength_vol*strength_vel, kind=wp) - kcomp(5)%sf(i, j, k)
740 t_kahan = updatedvar(5)%sf(i, j, k) + y_kahan
741 kcomp(5)%sf(i, j, k) = (t_kahan - updatedvar(5)%sf(i, j, k)) - y_kahan
742 updatedvar(5)%sf(i, j, k) = t_kahan
752# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
753#if defined(MFC_OpenACC)
754# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
756# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
757#elif defined(MFC_OpenMP)
758# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
760# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
762# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
771# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
773# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
775# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
777# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
779# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
781# 288 "/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"
806 real(wp),
dimension(3),
intent(in) :: center
807 integer,
dimension(3),
intent(in) :: cellaux
808 real(wp),
dimension(3),
intent(in) :: nodecoord
809 real(wp),
intent(in) :: stddsv
810 real(wp),
intent(in) :: strength_idx
811 real(wp),
intent(out) :: func
814 real(wp) :: theta, dtheta, L2, dzp, Lz2, zc
815 real(wp) :: Nr, Nr_count
817 distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + (center(3) - nodecoord(3))**2._wp)
819 if (num_dims == 3)
then
821 func = exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
827 nr = ceiling(2._wp*pi*nodecoord(2)/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
829 l2 = center(2)**2._wp + nodecoord(2)**2._wp - 2._wp*center(2)*nodecoord(2)*cos(theta)
830 distance = sqrt((center(1) - nodecoord(1))**2._wp + l2)
832 func = dtheta/2._wp/pi*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
834 do while (nr_count < nr - 1._wp)
835 nr_count = nr_count + 1._wp
836 theta = nr_count*dtheta
838 l2 = center(2)**2._wp + nodecoord(2)**2._wp - 2._wp*center(2)*nodecoord(2)*cos(theta)
839 distance = sqrt((center(1) - nodecoord(1))**2._wp + l2)
841 func = func + dtheta/2._wp/pi*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv) &
842 & **(3._wp*(strength_idx + 1._wp))
847 dzp = (lag_params%charwidth/(lag_params%charNz + 1._wp))
850 do i = 0, lag_params%charNz
851 zc = (-lag_params%charwidth/2._wp + dzp*(0.5_wp + i))
852 lz2 = (center(3) - zc)**2._wp
853 distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + lz2)
854 func = func + dzp/lag_params%charwidth*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
1353# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1355# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1357# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1359# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1361# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1363# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
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"
1388 real(wp),
intent(in) :: pos
1389 integer,
dimension(3),
intent(in) :: cell
1390 integer,
intent(in) :: i
1391 type(scalar_field),
dimension(sys_size),
intent(in) :: q_prim_vf
1393 real(wp),
dimension(5) :: xi, eta,
l
1395 if (fd_order == 2)
then
1397 xi(1) = x_cc(cell(1) - 1)
1398 eta(1) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 1, cell(2), cell(3))
1399 xi(2) = x_cc(cell(1))
1400 eta(2) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1), cell(2), cell(3))
1401 xi(3) = x_cc(cell(1) + 1)
1402 eta(3) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 1, cell(2), cell(3))
1403 else if (i == 2)
then
1404 xi(1) = y_cc(cell(2) - 1)
1405 eta(1) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 1, cell(3))
1406 xi(2) = y_cc(cell(2))
1407 eta(2) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2), cell(3))
1408 xi(3) = y_cc(cell(2) + 1)
1409 eta(3) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 1, cell(3))
1410 else if (i == 3)
then
1411 xi(1) = z_cc(cell(3) - 1)
1412 eta(1) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 1)
1413 xi(2) = z_cc(cell(3))
1414 eta(2) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3))
1415 xi(3) = z_cc(cell(3) + 1)
1416 eta(3) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 1)
1419 l(1) = ((pos - xi(2))*(pos - xi(3)))/((xi(1) - xi(2))*(xi(1) - xi(3)))
1420 l(2) = ((pos - xi(1))*(pos - xi(3)))/((xi(2) - xi(1))*(xi(2) - xi(3)))
1421 l(3) = ((pos - xi(1))*(pos - xi(2)))/((xi(3) - xi(1))*(xi(3) - xi(2)))
1423 v =
l(1)*eta(1) +
l(2)*eta(2) +
l(3)*eta(3)
1424 else if (fd_order == 4)
then
1426 xi(1) = x_cc(cell(1) - 2)
1427 eta(1) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 2, cell(2), cell(3))
1428 xi(2) = x_cc(cell(1) - 1)
1429 eta(2) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 1, cell(2), cell(3))
1430 xi(3) = x_cc(cell(1))
1431 eta(3) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1), cell(2), cell(3))
1432 xi(4) = x_cc(cell(1) + 1)
1433 eta(4) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 1, cell(2), cell(3))
1434 xi(5) = x_cc(cell(1) + 2)
1435 eta(5) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 2, cell(2), cell(3))
1436 else if (i == 2)
then
1437 xi(1) = y_cc(cell(2) - 2)
1438 eta(1) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 2, cell(3))
1439 xi(2) = y_cc(cell(2) - 1)
1440 eta(2) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 1, cell(3))
1441 xi(3) = y_cc(cell(2))
1442 eta(3) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2), cell(3))
1443 xi(4) = y_cc(cell(2) + 1)
1444 eta(4) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 1, cell(3))
1445 xi(5) = y_cc(cell(2) + 2)
1446 eta(5) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 2, cell(3))
1447 else if (i == 3)
then
1448 xi(1) = z_cc(cell(3) - 2)
1449 eta(1) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 2)
1450 xi(2) = z_cc(cell(3) - 1)
1451 eta(2) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 1)
1452 xi(3) = z_cc(cell(3))
1453 eta(3) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3))
1454 xi(4) = z_cc(cell(3) + 1)
1455 eta(4) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 1)
1456 xi(5) = z_cc(cell(3) + 2)
1457 eta(5) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 2)
1460 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)) &
1462 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)) &
1464 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)) &
1466 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)) &
1468 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)) &
1471 v =
l(1)*eta(1) +
l(2)*eta(2) +
l(3)*eta(3) +
l(4)*eta(4) +
l(5)*eta(5)
1481 function f_get_bubble_force(pos, rad, rdot, vel, mg, mv, Re, rho, cell, i, q_prim_vf)
result(force)
1484# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1486# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1488# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1490# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1492# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1494# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1496# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1498 real(wp),
intent(in) :: pos, rad, rdot, mg, mv, re, rho, vel
1499 integer,
dimension(3),
intent(in) :: cell
1500 integer,
intent(in) :: i
1501 type(scalar_field),
dimension(sys_size),
intent(in) :: q_prim_vf
1502 real(wp) :: dp, vol, force
1505 if (fd_order > 1)
then
1508 v_rel = vel - q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(cell(1), cell(2), cell(3))
1513 if (lag_params%drag_model == 1)
then
1514 force = force - (4._wp*pi*rad*v_rel)/re
1515 else if (lag_params%drag_model == 2)
then
1516 force = force - (6._wp*pi*rad*v_rel)/re
1517 else if (lag_params%drag_model == 3)
then
1518 force = force - (12._wp*pi*rad*v_rel)/re
1521 if (lag_pressure_force)
then
1524 dp =
grad_p_x(cell(1), cell(2), cell(3))
1525 else if (i == 2)
then
1526 dp =
grad_p_y(cell(1), cell(2), cell(3))
1527 else if (i == 3)
then
1528 dp =
grad_p_z(cell(1), cell(2), cell(3))
1531 vol = (4._wp/3._wp)*pi*(rad**3._wp)
1532 force = force - vol*dp
1535 if (lag_params%gravity_force)
then
1536 force = force + (mg + mv)*accel_bf(i)