Skip to content

Commit 96d5f2f

Browse files
wilfonbawilfonba
committed
remove IB support and precompute pressure gradients
1 parent b0939fd commit 96d5f2f

2 files changed

Lines changed: 141 additions & 94 deletions

File tree

‎src/simulation/m_bubbles_EL.fpp‎

Lines changed: 54 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,8 @@ module m_bubbles_EL
3131

3232
use m_ibm
3333

34+
use m_finite_differences
35+
3436
implicit none
3537

3638
!(nBub)
@@ -189,6 +191,29 @@ contains
189191
$:GPU_UPDATE(device='[moving_lag_bubbles, lag_pressure_force, &
190192
& lag_gravity_force, lag_vel_model, lag_drag_model]')
191193
194+
! Allocate cell-centered pressure gradient arrays and FD coefficients
195+
if (lag_params%vel_model > 0 .and. lag_params%pressure_force) then
196+
@:ALLOCATE(grad_p_x(0:m, 0:n, 0:p))
197+
@:ALLOCATE(fd_coeff_x_pgrad(-fd_number:fd_number, 0:m))
198+
call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_pgrad, &
199+
buff_size, fd_number, fd_order)
200+
$:GPU_UPDATE(device='[fd_coeff_x_pgrad]')
201+
if (n > 0) then
202+
@:ALLOCATE(grad_p_y(0:m, 0:n, 0:p))
203+
@:ALLOCATE(fd_coeff_y_pgrad(-fd_number:fd_number, 0:n))
204+
call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_pgrad, &
205+
buff_size, fd_number, fd_order)
206+
$:GPU_UPDATE(device='[fd_coeff_y_pgrad]')
207+
end if
208+
if (p > 0) then
209+
@:ALLOCATE(grad_p_z(0:m, 0:n, 0:p))
210+
@:ALLOCATE(fd_coeff_z_pgrad(-fd_number:fd_number, 0:p))
211+
call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_pgrad, &
212+
buff_size, fd_number, fd_order)
213+
$:GPU_UPDATE(device='[fd_coeff_z_pgrad]')
214+
end if
215+
end if
216+
192217
call s_read_input_bubbles(q_cons_vf, bc_type)
193218
194219
end subroutine s_initialize_bubbles_EL_module
@@ -612,10 +637,9 @@ contains
612637

613638
integer :: k, l
614639

615-
call nvtxStartRange("LAGRANGE-BUBBLE-DYNAMICS")
616-
617640
! Subgrid p_inf model based on Maeda and Colonius (2018).
618641
if (lag_params%pressure_corrector) then
642+
call nvtxStartRange("LAGRANGE-BUBBLE-PINF-CORRECTION")
619643
! Calculate velocity potentials (valid for one bubble per cell)
620644
$:GPU_PARALLEL_LOOP(private='[k,cell,paux,preterm1,term2,Romega,myR0,myR,myV,myPb,pint,term1_fac]')
621645
do k = 1, n_el_bubs_loc
@@ -635,8 +659,17 @@ contains
635659
end if
636660
end do
637661
$:END_GPU_PARALLEL_LOOP()
662+
call nvtxEndRange()
663+
end if
664+
665+
! Precompute cell-centered pressure gradients for translational motion
666+
if (moving_lag_bubbles .and. lag_pressure_force) then
667+
call nvtxStartRange("LAGRANGE-BUBBLE-PRESSURE-GRADIENT")
668+
call s_compute_pressure_gradients(q_prim_vf)
669+
call nvtxEndRange()
638670
end if
639671

672+
call nvtxStartRange("LAGRANGE-BUBBLE-DYNAMICS")
640673
! Radial motion model
641674
adap_dt_stop_sum = 0
642675
$:GPU_PARALLEL_LOOP(private='[k,myalpha_rho,myalpha,Re,cell,myVapFlux,preterm1, term2, paux, pint, Romega,term1_fac,myR_m, mygamma_m, myPb, myMass_n, myMass_v,myR, myV, myBeta_c, myBeta_t, myR0, myPbdot, myMvdot,myPinf, aux1,aux2, myCson, myRho,gamma,pi_inf,qv,dmalf, dmntait, dmBtait, dm_bub_adv_src, dm_divu,adap_dt_stop,myPos,myVel]', &
@@ -737,6 +770,7 @@ contains
737770

738771
end do
739772
$:END_GPU_PARALLEL_LOOP()
773+
call nvtxEndRange
740774

741775
if (adap_dt .and. adap_dt_stop_sum > 0) call s_mpi_abort("Adaptive time stepping failed to converge.")
742776

@@ -746,8 +780,6 @@ contains
746780
call s_smear_voidfraction(bc_type)
747781
end if
748782

749-
call nvtxEndRange
750-
751783
end subroutine s_compute_bubble_EL_dynamics
752784

753785
!> The purpose of this subroutine is to obtain the bubble source terms based on Maeda and Colonius (2018)
@@ -763,6 +795,7 @@ contains
763795

764796
integer :: i, j, k, l
765797

798+
call nvtxStartRange("LAGRANGE-BUBBLE-EL-SOURCE")
766799
! (q / (1 - beta)) * d(beta)/dt source
767800
if (p == 0) then
768801
$:GPU_PARALLEL_LOOP(private='[i,j,k,l]', collapse=4)
@@ -846,6 +879,7 @@ contains
846879
end do
847880
$:END_GPU_PARALLEL_LOOP()
848881
end do
882+
call nvtxEndRange ! LAGRANGE-BUBBLE-EL-SOURCE
849883

850884
end subroutine s_compute_bubbles_EL_source
851885

@@ -892,7 +926,7 @@ contains
892926
type(integer_field), dimension(1:num_dims, 1:2), intent(in) :: bc_type
893927
integer :: i, j, k, l
894928

895-
call nvtxStartRange("BUBBLES-LAGRANGE-KERNELS")
929+
call nvtxStartRange("BUBBLES-LAGRANGE-SMEARING")
896930
$:GPU_PARALLEL_LOOP(private='[i,j,k,l]', collapse=4)
897931
do i = 1, q_beta_idx
898932
do l = idwbuff(3)%beg, idwbuff(3)%end
@@ -929,7 +963,7 @@ contains
929963
end do
930964
end do
931965
$:END_GPU_PARALLEL_LOOP()
932-
call nvtxEndRange ! BUBBLES-LAGRANGE-KERNELS
966+
call nvtxEndRange ! BUBBLES-LAGRANGE-SMEARING
933967

934968
end subroutine s_smear_voidfraction
935969

@@ -1483,26 +1517,6 @@ contains
14831517
if (q_prim_vf(advxb)%sf(cell(1), cell(2), cell(3)) < (1._wp - lag_params%valmaxvoid)) then
14841518
keep_bubble(k) = 0
14851519
end if
1486-
1487-
! Move bubbles back to surface of IB
1488-
if (ib) then
1489-
cell = fd_number - buff_size
1490-
call s_locate_cell(mtn_pos(k, 1:3, 2), cell, mtn_s(k, 1:3, 2))
1491-
1492-
if (ib_markers%sf(cell(1), cell(2), cell(3)) /= 0) then
1493-
patch_id = ib_markers%sf(cell(1), cell(2), cell(3))
1494-
1495-
$:GPU_LOOP(parallelism='[seq]')
1496-
do i = 1, num_dims
1497-
mtn_pos(k, i, 2) = mtn_pos(k, i, 2) - &
1498-
levelset_norm%sf(cell(1), cell(2), cell(3), patch_id, i) &
1499-
*levelset%sf(cell(1), cell(2), cell(3), patch_id)
1500-
end do
1501-
1502-
cell = fd_number - buff_size
1503-
call s_locate_cell(mtn_pos(k, 1:3, 2), cell, mtn_s(k, 1:3, 2))
1504-
end if
1505-
end if
15061520
end if
15071521
end do
15081522
$:END_GPU_PARALLEL_LOOP()
@@ -2229,6 +2243,20 @@ contains
22292243
@:DEALLOCATE(mtn_dposdt)
22302244
@:DEALLOCATE(mtn_dveldt)
22312245

2246+
! Deallocate pressure gradient arrays and FD coefficients
2247+
if (lag_params%vel_model > 0 .and. lag_params%pressure_force) then
2248+
@:DEALLOCATE(grad_p_x)
2249+
@:DEALLOCATE(fd_coeff_x_pgrad)
2250+
if (n > 0) then
2251+
@:DEALLOCATE(grad_p_y)
2252+
@:DEALLOCATE(fd_coeff_y_pgrad)
2253+
if (p > 0) then
2254+
@:DEALLOCATE(grad_p_z)
2255+
@:DEALLOCATE(fd_coeff_z_pgrad)
2256+
end if
2257+
end if
2258+
end if
2259+
22322260
end subroutine s_finalize_lagrangian_solver
22332261

22342262
end module m_bubbles_EL

‎src/simulation/m_bubbles_EL_kernels.fpp‎

Lines changed: 87 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,16 @@ module m_bubbles_EL_kernels
1212

1313
implicit none
1414

15+
! Cell-centered pressure gradients (precomputed for translational motion)
16+
real(wp), allocatable, dimension(:, :, :) :: grad_p_x, grad_p_y, grad_p_z
17+
$:GPU_DECLARE(create='[grad_p_x, grad_p_y, grad_p_z]')
18+
19+
! Finite-difference coefficients for pressure gradient computation
20+
real(wp), allocatable, dimension(:, :) :: fd_coeff_x_pgrad
21+
real(wp), allocatable, dimension(:, :) :: fd_coeff_y_pgrad
22+
real(wp), allocatable, dimension(:, :) :: fd_coeff_z_pgrad
23+
$:GPU_DECLARE(create='[fd_coeff_x_pgrad, fd_coeff_y_pgrad, fd_coeff_z_pgrad]')
24+
1525
contains
1626

1727
!> The purpose of this subroutine is to smear the strength of the lagrangian
@@ -372,6 +382,71 @@ contains
372382

373383
end subroutine s_get_cell
374384

385+
!> Precomputes cell-centered pressure gradients (dp/dx, dp/dy, dp/dz) at all cell centers
386+
!! using finite-difference coefficients of the specified order. This avoids
387+
!! scattered memory accesses to the pressure field when computing translational
388+
!! bubble forces.
389+
!! @param q_prim_vf Primitive variables (pressure is at index E_idx)
390+
subroutine s_compute_pressure_gradients(q_prim_vf)
391+
392+
type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
393+
394+
integer :: i, j, k, r
395+
396+
! dp/dx at all cell centers
397+
$:GPU_PARALLEL_LOOP(private='[i,j,k,r]', collapse=3)
398+
do k = 0, p
399+
do j = 0, n
400+
do i = 0, m
401+
grad_p_x(i, j, k) = 0._wp
402+
$:GPU_LOOP(parallelism='[seq]')
403+
do r = -fd_number, fd_number
404+
grad_p_x(i, j, k) = grad_p_x(i, j, k) + &
405+
q_prim_vf(E_idx)%sf(i + r, j, k)*fd_coeff_x_pgrad(r, i)
406+
end do
407+
end do
408+
end do
409+
end do
410+
$:END_GPU_PARALLEL_LOOP()
411+
412+
! dp/dy at all cell centers
413+
if (n > 0) then
414+
$:GPU_PARALLEL_LOOP(private='[i,j,k,r]', collapse=3)
415+
do k = 0, p
416+
do j = 0, n
417+
do i = 0, m
418+
grad_p_y(i, j, k) = 0._wp
419+
$:GPU_LOOP(parallelism='[seq]')
420+
do r = -fd_number, fd_number
421+
grad_p_y(i, j, k) = grad_p_y(i, j, k) + &
422+
q_prim_vf(E_idx)%sf(i, j + r, k)*fd_coeff_y_pgrad(r, j)
423+
end do
424+
end do
425+
end do
426+
end do
427+
$:END_GPU_PARALLEL_LOOP()
428+
end if
429+
430+
! dp/dz at all cell centers
431+
if (p > 0) then
432+
$:GPU_PARALLEL_LOOP(private='[i,j,k,r]', collapse=3)
433+
do k = 0, p
434+
do j = 0, n
435+
do i = 0, m
436+
grad_p_z(i, j, k) = 0._wp
437+
$:GPU_LOOP(parallelism='[seq]')
438+
do r = -fd_number, fd_number
439+
grad_p_z(i, j, k) = grad_p_z(i, j, k) + &
440+
q_prim_vf(E_idx)%sf(i, j, k + r)*fd_coeff_z_pgrad(r, k)
441+
end do
442+
end do
443+
end do
444+
end do
445+
$:END_GPU_PARALLEL_LOOP()
446+
end if
447+
448+
end subroutine s_compute_pressure_gradients
449+
375450
!! This function interpolates the velocity of Eulerian field at the position
376451
!! of the bubble.
377452
!! @param pos Position of the bubble in directiion i
@@ -491,74 +566,8 @@ contains
491566
integer, intent(in) :: i
492567
type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
493568

494-
real(wp) :: a, dp, vol, force
569+
real(wp) :: dp, vol, force
495570
real(wp) :: v_rel
496-
real(wp), dimension(fd_order) :: xi, eta, L
497-
498-
if (fd_order <= 2) then
499-
if (i == 1) then
500-
dp = (q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3)) - &
501-
q_prim_vf(E_idx)%sf(cell(1) - 1, cell(2), cell(3)))/ &
502-
(x_cc(cell(1) + 1) - x_cc(cell(1) - 1))
503-
elseif (i == 2) then
504-
dp = (q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3)) - &
505-
q_prim_vf(E_idx)%sf(cell(1), cell(2) - 1, cell(3)))/ &
506-
(y_cc(cell(2) + 1) - y_cc(cell(2) - 1))
507-
elseif (i == 3) then
508-
dp = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1) - &
509-
q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) - 1))/ &
510-
(z_cc(cell(3) + 1) - z_cc(cell(3) - 1))
511-
end if
512-
elseif (fd_order == 4) then
513-
if (i == 1) then
514-
xi(1) = x_cc(cell(1) - 1)
515-
eta(1) = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)) - &
516-
q_prim_vf(E_idx)%sf(cell(1) - 2, cell(2), cell(3)))/ &
517-
(x_cc(cell(1)) - x_cc(cell(1) - 2))
518-
xi(2) = x_cc(cell(1))
519-
eta(2) = (q_prim_vf(E_idx)%sf(cell(1) + 1, cell(2), cell(3)) - &
520-
q_prim_vf(E_idx)%sf(cell(1) - 1, cell(2), cell(3)))/ &
521-
(x_cc(cell(1) + 1) - x_cc(cell(1) - 1))
522-
xi(3) = x_cc(cell(1) + 1)
523-
eta(3) = (q_prim_vf(E_idx)%sf(cell(1) + 2, cell(2), cell(3)) - &
524-
q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)))/ &
525-
(x_cc(cell(1) + 2) - x_cc(cell(1)))
526-
elseif (i == 2) then
527-
xi(1) = y_cc(cell(2) - 1)
528-
eta(1) = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)) - &
529-
q_prim_vf(E_idx)%sf(cell(1), cell(2) - 2, cell(3)))/ &
530-
(y_cc(cell(2)) - y_cc(cell(2) - 2))
531-
xi(2) = y_cc(cell(2))
532-
eta(2) = (q_prim_vf(E_idx)%sf(cell(1), cell(2) + 1, cell(3)) - &
533-
q_prim_vf(E_idx)%sf(cell(1), cell(2) - 1, cell(3)))/ &
534-
(y_cc(cell(2) + 1) - y_cc(cell(2) - 1))
535-
xi(3) = y_cc(cell(2) + 1)
536-
eta(3) = (q_prim_vf(E_idx)%sf(cell(1), cell(2) + 2, cell(3)) - &
537-
q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)))/ &
538-
(y_cc(cell(2) + 2) - y_cc(cell(2)))
539-
elseif (i == 3) then
540-
xi(1) = z_cc(cell(3) - 1)
541-
eta(1) = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)) - &
542-
q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) - 2))/ &
543-
(z_cc(cell(3)) - z_cc(cell(3) - 2))
544-
xi(2) = z_cc(cell(3))
545-
eta(2) = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 1) - &
546-
q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) - 1))/ &
547-
(z_cc(cell(3) + 1) - z_cc(cell(3) - 1))
548-
xi(3) = z_cc(cell(3) + 1)
549-
eta(3) = (q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3) + 2) - &
550-
q_prim_vf(E_idx)%sf(cell(1), cell(2), cell(3)))/ &
551-
(z_cc(cell(3) + 2) - z_cc(cell(3)))
552-
end if
553-
554-
L(1) = ((pos - xi(2))*(pos - xi(3)))/((xi(1) - xi(2))*(xi(1) - xi(3)))
555-
L(2) = ((pos - xi(1))*(pos - xi(3)))/((xi(2) - xi(1))*(xi(2) - xi(3)))
556-
L(3) = ((pos - xi(1))*(pos - xi(2)))/((xi(3) - xi(1))*(xi(3) - xi(2)))
557-
558-
dp = L(1)*eta(1) + L(2)*eta(2) + L(3)*eta(3)
559-
end if
560-
561-
vol = (4._wp/3._wp)*pi*(rad**3._wp)
562571

563572
if (fd_order > 1) then
564573
v_rel = vel - f_interpolate_velocity(pos, cell, i, q_prim_vf)
@@ -576,7 +585,17 @@ contains
576585
force = force - (12._wp*pi*rad*v_rel)/Re
577586
end if
578587

579-
if (lag_params%pressure_force) then
588+
if (lag_pressure_force) then
589+
! Use precomputed cell-centered pressure gradients
590+
if (i == 1) then
591+
dp = grad_p_x(cell(1), cell(2), cell(3))
592+
elseif (i == 2) then
593+
dp = grad_p_y(cell(1), cell(2), cell(3))
594+
elseif (i == 3) then
595+
dp = grad_p_z(cell(1), cell(2), cell(3))
596+
end if
597+
598+
vol = (4._wp/3._wp)*pi*(rad**3._wp)
580599
force = force - vol*dp
581600
end if
582601

0 commit comments

Comments
 (0)