Skip to content

Commit 45d98f0

Browse files
danieljvickersdanieljvickers
committed
Initial commit of higher-order central differencing for the IB force volume integrals
1 parent 8107b97 commit 45d98f0

3 files changed

Lines changed: 39 additions & 65 deletions

File tree

‎src/simulation/m_derived_variables.fpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ contains
5050
! to be implemented in the subroutine s_compute_finite_difference_coefficients.
5151

5252
! Allocating centered finite-difference coefficients
53-
if (probe_wrt) then
53+
if (probe_wrt .or. ib) then
5454
@:ALLOCATE(fd_coeff_x(-fd_number:fd_number, 0:m))
5555
if (n > 0) then
5656
@:ALLOCATE(fd_coeff_y(-fd_number:fd_number, 0:n))
@@ -74,9 +74,9 @@ contains
7474
!> Allocate and open derived variables. Computing FD coefficients.
7575
impure subroutine s_initialize_derived_variables
7676

77-
if (probe_wrt) then
77+
if (probe_wrt .or. ib) then
7878
! Opening and writing header of flow probe files
79-
if (proc_rank == 0) then
79+
if (proc_rank == 0 .and. probe_wrt) then
8080
call s_open_probe_files()
8181
call s_open_com_files()
8282
end if

‎src/simulation/m_global_parameters.fpp‎

Lines changed: 1 addition & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -818,15 +818,7 @@ contains
818818

819819
if (ib) allocate (MPI_IO_IB_DATA%var%sf(0:m,0:n,0:p))
820820

821-
if (elasticity) then
822-
fd_number = max(1, fd_order/2)
823-
end if
824-
825-
if (mhd) then
826-
fd_number = max(1, fd_order/2)
827-
end if
828-
829-
if (probe_wrt) then
821+
if (elasticity .or. mhd .or. probe_wrt .or. ib) then
830822
fd_number = max(1, fd_order/2)
831823
end if
832824

‎src/simulation/m_ibm.fpp‎

Lines changed: 35 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
module m_ibm
1010

1111
use m_derived_types
12+
use m_derived_variables
1213
use m_global_parameters
1314
use m_mpi_proxy
1415
use m_variables_conversion
@@ -910,7 +911,7 @@ contains
910911
integer :: i, j, k, l, encoded_ib_idx, ib_idx, ib_idx_temp, fluid_idx
911912
real(wp), dimension(num_ibs, 3) :: forces, torques
912913
! viscous stress tensor with temp vectors to hold divergence calculations
913-
real(wp), dimension(1:3,1:3) :: viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2
914+
real(wp), dimension(1:3,1:3) :: viscous_stress
914915
real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution
915916
real(wp) :: cell_volume, dx, dy, dz, dynamic_viscosity
916917

@@ -936,9 +937,8 @@ contains
936937
end if
937938

938939
$:GPU_PARALLEL_LOOP(private='[ib_idx, ib_idx_temp, encoded_ib_idx, fluid_idx, radial_vector, local_force_contribution, &
939-
& cell_volume, local_torque_contribution, dynamic_viscosity, viscous_stress_div, &
940-
& viscous_stress_div_1, viscous_stress_div_2, dx, dy, dz]', copy='[forces, torques]', &
941-
& copyin='[dynamic_viscosities]', collapse=3)
940+
& cell_volume, local_torque_contribution, dynamic_viscosity, viscous_stress]', copy='[forces, &
941+
& torques]', copyin='[dynamic_viscosities]', collapse=3)
942942
do i = 0, m
943943
do j = 0, n
944944
do k = 0, p
@@ -955,28 +955,21 @@ contains
955955
radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, &
956956
& patch_ib(ib_idx)%y_centroid, 0._wp]
957957
end if
958-
dx = x_cc(i + 1) - x_cc(i)
959-
dy = y_cc(j + 1) - y_cc(j)
960958

961959
local_force_contribution(:) = 0._wp
960+
961+
! compute the pressure force component, which is the negative pressure gradient
962962
do fluid_idx = 0, num_fluids - 1
963-
! Get the pressure contribution to force via a finite difference to compute the 2D components of the
964-
! gradient of the pressure and cell volume
965-
local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(eqn_idx%E &
966-
& + fluid_idx)%sf(i + 1, j, &
967-
& k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i - 1, j, k))/(2._wp*dx) ! force is the negative pressure gradient
968-
local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(eqn_idx%E &
969-
& + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, &
970-
& j - 1, k))/(2._wp*dy)
971-
cell_volume = abs(dx*dy)
972-
! add the 3D component of the pressure gradient, if we are working in 3 dimensions
973-
if (num_dims == 3) then
974-
dz = z_cc(k + 1) - z_cc(k)
975-
local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(eqn_idx%E &
976-
& + fluid_idx)%sf(i, j, &
977-
& k + 1) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j, k - 1))/(2._wp*dz)
978-
cell_volume = abs(cell_volume*dz)
979-
end if
963+
do l = -fd_number, fd_number
964+
local_force_contribution(1) = local_force_contribution(1) - (fd%fd_coeff_x(l, &
965+
& i)*q_prim_vf(eqn_idx%E + fluid_idx)%sf(i + l, j, k))
966+
local_force_contribution(2) = local_force_contribution(2) - (fd%fd_coeff_y(l, &
967+
& j)*q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j + l, k))
968+
if (num_dims == 3) then
969+
local_force_contribution(3) = local_force_contribution(3) - (fd%fd_coeff_z(l, &
970+
& k)*q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j, k + l))
971+
end if
972+
end do
980973
end do
981974

982975
! get the viscous stress and add its contribution if that is considered
@@ -989,41 +982,30 @@ contains
989982
& k)*dynamic_viscosities(fluid_idx))
990983
end do
991984

992-
! get the linear force components first
993-
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, &
994-
& j, k)
995-
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, &
996-
& j, k)
997-
! get x derivative of the first-row of viscous stress tensor
998-
viscous_stress_div(1,1:3) = (viscous_stress_div_2(1,1:3) - viscous_stress_div_1(1,1:3))/(2._wp*dx)
999-
! add the x components of the divergence to the force
1000-
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1,1:3)
1001-
1002-
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, &
1003-
& j - 1, k)
1004-
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, &
1005-
& j + 1, k)
1006-
! get y derivative of the second-row of viscous stress tensor
1007-
viscous_stress_div(2,1:3) = (viscous_stress_div_2(2,1:3) - viscous_stress_div_1(2,1:3))/(2._wp*dy)
1008-
! add the y components of the divergence to the force
1009-
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2,1:3)
1010-
1011-
if (num_dims == 3) then
1012-
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, &
1013-
& j, k - 1)
1014-
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, &
1015-
& j, k + 1)
1016-
viscous_stress_div(3,1:3) = (viscous_stress_div_2(3,1:3) - viscous_stress_div_1(3, &
1017-
& 1:3))/(2._wp*dz)
1018-
! add the z components of the divergence to the force
1019-
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3,1:3)
1020-
end if
985+
do l = -fd_number, fd_number
986+
call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i + l, j, k)
987+
local_force_contribution(1:3) = local_force_contribution(1:3) + fd%fd_coeff_x(l, &
988+
& i)*viscous_stress(1,1:3)
989+
990+
call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i, j + l, k)
991+
local_force_contribution(1:3) = local_force_contribution(1:3) + fd%fd_coeff_x(l, &
992+
& i)*viscous_stress(2,1:3)
993+
994+
if (num_dims == 3) then
995+
call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i, j, &
996+
& k + l)
997+
local_force_contribution(1:3) = local_force_contribution(1:3) + fd%fd_coeff_x(l, &
998+
& i)*viscous_stress(3,1:3)
999+
end if
1000+
end do
10211001
end if
10221002

10231003
call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution)
10241004

10251005
! Update the force and torque values atomically to prevent race conditions
1026-
do l = 1, 3
1006+
cell_volume = dx(i)*dy(j)
1007+
if (num_dims == 3) cell_volume = cell_volume*dz(k)
1008+
do l = 1, num_dims
10271009
$:GPU_ATOMIC(atomic='update')
10281010
forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume)
10291011
$:GPU_ATOMIC(atomic='update')

0 commit comments

Comments
 (0)