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"
625# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
626#elif defined(MFC_OpenMP)
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# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
641# 210 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
642 do k = smear_z_beg, smear_z_end
643 do j = smear_y_beg, smear_y_end
644 do i = smear_x_beg, smear_x_end
649 nodecoord(1) = x_cc(i)
650 nodecoord(2) = y_cc(j)
652 if (p > 0) nodecoord(3) = z_cc(k)
655 di_beg = max(i - mapcells, 0)
656 di_end = min(i + mapcells, m)
657 dj_beg = max(j - mapcells, 0)
658 dj_end = min(j + mapcells, n)
659 dk_beg = max(k - mapcells, 0)
660 dk_end = min(k + mapcells, p)
663# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
664#if defined(MFC_OpenACC)
665# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
667# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
668#elif defined(MFC_OpenMP)
669# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
671# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
673 do dk = dk_beg, dk_end
675# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
676#if defined(MFC_OpenACC)
677# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
679# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
680#elif defined(MFC_OpenMP)
681# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
683# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
685 do dj = dj_beg, dj_end
687# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
688#if defined(MFC_OpenACC)
689# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
691# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
692#elif defined(MFC_OpenMP)
693# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
695# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
697 do di = di_beg, di_end
699# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
700#if defined(MFC_OpenACC)
701# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
703# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
704#elif defined(MFC_OpenMP)
705# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
707# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
713 volpart = 4._wp/3._wp*pi*lbk_rad(bub_idx, 2)**3._wp
714 s_coord(1:3) = lbk_s(bub_idx,1:3,2)
718 strength_vol = volpart
719 strength_vel = 4._wp*pi*lbk_rad(bub_idx, 2)**2._wp*lbk_vel(bub_idx, 2)
721 center(1:2) = lbk_pos(bub_idx,1:2,2)
723 if (p > 0) center(3) = lbk_pos(bub_idx, 3, 2)
728 y_kahan = real(func*strength_vol, kind=wp) - kcomp(1)%sf(i, j, k)
729 t_kahan = updatedvar(1)%sf(i, j, k) + y_kahan
730 kcomp(1)%sf(i, j, k) = (t_kahan - updatedvar(1)%sf(i, j, k)) - y_kahan
731 updatedvar(1)%sf(i, j, k) = t_kahan
734 y_kahan = real(func*strength_vel, kind=wp) - kcomp(2)%sf(i, j, k)
735 t_kahan = updatedvar(2)%sf(i, j, k) + y_kahan
736 kcomp(2)%sf(i, j, k) = (t_kahan - updatedvar(2)%sf(i, j, k)) - y_kahan
737 updatedvar(2)%sf(i, j, k) = t_kahan
739 if (lag_params%cluster_type >= 4)
then
741 y_kahan = real(func2*strength_vol*strength_vel, kind=wp) - kcomp(5)%sf(i, j, k)
742 t_kahan = updatedvar(5)%sf(i, j, k) + y_kahan
743 kcomp(5)%sf(i, j, k) = (t_kahan - updatedvar(5)%sf(i, j, k)) - y_kahan
744 updatedvar(5)%sf(i, j, k) = t_kahan
754# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
755#if defined(MFC_OpenACC)
756# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
758# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
759#elif defined(MFC_OpenMP)
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"
764# 281 "/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"
805# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
808 real(wp),
dimension(3),
intent(in) :: center
809 integer,
dimension(3),
intent(in) :: cellaux
810 real(wp),
dimension(3),
intent(in) :: nodecoord
811 real(wp),
intent(in) :: stddsv
812 real(wp),
intent(in) :: strength_idx
813 real(wp),
intent(out) :: func
816 real(wp) :: theta, dtheta, L2, dzp, Lz2, zc
817 real(wp) :: Nr, Nr_count
819 distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + (center(3) - nodecoord(3))**2._wp)
821 if (num_dims == 3)
then
823 func = exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
829 nr = ceiling(2._wp*pi*nodecoord(2)/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
831 l2 = center(2)**2._wp + nodecoord(2)**2._wp - 2._wp*center(2)*nodecoord(2)*cos(theta)
832 distance = sqrt((center(1) - nodecoord(1))**2._wp + l2)
834 func = dtheta/2._wp/pi*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
836 do while (nr_count < nr - 1._wp)
837 nr_count = nr_count + 1._wp
838 theta = nr_count*dtheta
840 l2 = center(2)**2._wp + nodecoord(2)**2._wp - 2._wp*center(2)*nodecoord(2)*cos(theta)
841 distance = sqrt((center(1) - nodecoord(1))**2._wp + l2)
843 func = func + dtheta/2._wp/pi*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv) &
844 & **(3._wp*(strength_idx + 1._wp))
849 dzp = (lag_params%charwidth/(lag_params%charNz + 1._wp))
852 do i = 0, lag_params%charNz
853 zc = (-lag_params%charwidth/2._wp + dzp*(0.5_wp + i))
854 lz2 = (center(3) - zc)**2._wp
855 distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + lz2)
856 func = func + dzp/lag_params%charwidth*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
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"
1387# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1390 real(wp),
intent(in) :: pos
1391 integer,
dimension(3),
intent(in) :: cell
1392 integer,
intent(in) :: i
1393 type(scalar_field),
dimension(sys_size),
intent(in) :: q_prim_vf
1395 real(wp),
dimension(5) :: xi, eta,
l
1397 if (fd_order == 2)
then
1399 xi(1) = x_cc(cell(1) - 1)
1400 eta(1) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 1, cell(2), cell(3))
1401 xi(2) = x_cc(cell(1))
1402 eta(2) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1), cell(2), cell(3))
1403 xi(3) = x_cc(cell(1) + 1)
1404 eta(3) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 1, cell(2), cell(3))
1405 else if (i == 2)
then
1406 xi(1) = y_cc(cell(2) - 1)
1407 eta(1) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 1, cell(3))
1408 xi(2) = y_cc(cell(2))
1409 eta(2) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2), cell(3))
1410 xi(3) = y_cc(cell(2) + 1)
1411 eta(3) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 1, cell(3))
1412 else if (i == 3)
then
1413 xi(1) = z_cc(cell(3) - 1)
1414 eta(1) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 1)
1415 xi(2) = z_cc(cell(3))
1416 eta(2) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3))
1417 xi(3) = z_cc(cell(3) + 1)
1418 eta(3) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 1)
1421 l(1) = ((pos - xi(2))*(pos - xi(3)))/((xi(1) - xi(2))*(xi(1) - xi(3)))
1422 l(2) = ((pos - xi(1))*(pos - xi(3)))/((xi(2) - xi(1))*(xi(2) - xi(3)))
1423 l(3) = ((pos - xi(1))*(pos - xi(2)))/((xi(3) - xi(1))*(xi(3) - xi(2)))
1425 v =
l(1)*eta(1) +
l(2)*eta(2) +
l(3)*eta(3)
1426 else if (fd_order == 4)
then
1428 xi(1) = x_cc(cell(1) - 2)
1429 eta(1) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 2, cell(2), cell(3))
1430 xi(2) = x_cc(cell(1) - 1)
1431 eta(2) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 1, cell(2), cell(3))
1432 xi(3) = x_cc(cell(1))
1433 eta(3) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1), cell(2), cell(3))
1434 xi(4) = x_cc(cell(1) + 1)
1435 eta(4) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 1, cell(2), cell(3))
1436 xi(5) = x_cc(cell(1) + 2)
1437 eta(5) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 2, cell(2), cell(3))
1438 else if (i == 2)
then
1439 xi(1) = y_cc(cell(2) - 2)
1440 eta(1) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 2, cell(3))
1441 xi(2) = y_cc(cell(2) - 1)
1442 eta(2) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 1, cell(3))
1443 xi(3) = y_cc(cell(2))
1444 eta(3) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2), cell(3))
1445 xi(4) = y_cc(cell(2) + 1)
1446 eta(4) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 1, cell(3))
1447 xi(5) = y_cc(cell(2) + 2)
1448 eta(5) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 2, cell(3))
1449 else if (i == 3)
then
1450 xi(1) = z_cc(cell(3) - 2)
1451 eta(1) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 2)
1452 xi(2) = z_cc(cell(3) - 1)
1453 eta(2) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 1)
1454 xi(3) = z_cc(cell(3))
1455 eta(3) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3))
1456 xi(4) = z_cc(cell(3) + 1)
1457 eta(4) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 1)
1458 xi(5) = z_cc(cell(3) + 2)
1459 eta(5) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 2)
1462 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)) &
1464 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)) &
1466 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)) &
1468 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)) &
1470 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)) &
1473 v =
l(1)*eta(1) +
l(2)*eta(2) +
l(3)*eta(3) +
l(4)*eta(4) +
l(5)*eta(5)
1483 function f_get_bubble_force(pos, rad, rdot, vel, mg, mv, Re, rho, cell, i, q_prim_vf)
result(force)
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# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1500 real(wp),
intent(in) :: pos, rad, rdot, mg, mv, re, rho, vel
1501 integer,
dimension(3),
intent(in) :: cell
1502 integer,
intent(in) :: i
1503 type(scalar_field),
dimension(sys_size),
intent(in) :: q_prim_vf
1504 real(wp) :: dp, vol, force
1507 if (fd_order > 1)
then
1510 v_rel = vel - q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(cell(1), cell(2), cell(3))
1515 if (lag_params%drag_model == 1)
then
1516 force = force - (4._wp*pi*rad*v_rel)/re
1517 else if (lag_params%drag_model == 2)
then
1518 force = force - (6._wp*pi*rad*v_rel)/re
1519 else if (lag_params%drag_model == 3)
then
1520 force = force - (12._wp*pi*rad*v_rel)/re
1523 if (lag_pressure_force)
then
1526 dp =
grad_p_x(cell(1), cell(2), cell(3))
1527 else if (i == 2)
then
1528 dp =
grad_p_y(cell(1), cell(2), cell(3))
1529 else if (i == 3)
then
1530 dp =
grad_p_z(cell(1), cell(2), cell(3))
1533 vol = (4._wp/3._wp)*pi*(rad**3._wp)
1534 force = force - vol*dp
1537 if (lag_params%gravity_force)
then
1538 force = force + (mg + mv)*accel_bf(i)