Skip to content

Fix slip BCs for IBM #789

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Jan 19, 2025
Merged

Fix slip BCs for IBM #789

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b0b1d76
pull fixes from Henrys pr
sbryngelson Jan 15, 2025
a45c4dc
Update submit.sh
sbryngelson Jan 18, 2025
a9a0f3c
Update submit-bench.sh
sbryngelson Jan 18, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/common/m_constants.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ module m_constants
integer, parameter :: num_patches_max = 10
integer, parameter :: pathlen_max = 400
integer, parameter :: nnode = 4 !< Number of QBMM nodes
integer, parameter :: gp_layers = 3 !< Number of ghost point layers for IBM
real(wp), parameter :: capillary_cutoff = 1e-6 !< color function gradient magnitude at which to apply the surface tension fluxes
real(wp), parameter :: acoustic_spatial_support_width = 2.5_wp !< Spatial support width of acoustic source, used in s_source_spatial
real(wp), parameter :: dflt_vcfl_dt = 100._wp !< value of vcfl_dt when viscosity is off for computing adaptive timestep size
Expand Down
4 changes: 1 addition & 3 deletions src/common/m_derived_types.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -367,15 +367,13 @@ module m_derived_types

!> Ghost Point for Immersed Boundaries
type ghost_point

real(wp), dimension(3) :: loc !< Physical location of the ghost point
integer, dimension(3) :: loc !< Physical location of the ghost point
real(wp), dimension(3) :: ip_loc !< Physical location of the image point
integer, dimension(3) :: ip_grid !< Top left grid point of IP
real(wp), dimension(2, 2, 2) :: interp_coeffs !< Interpolation Coefficients of image point
integer :: ib_patch_id !< ID of the IB Patch the ghost point is part of
logical :: slip
integer, dimension(3) :: DB

end type ghost_point

!> Species parameters
Expand Down
13 changes: 7 additions & 6 deletions src/simulation/m_ibm.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@

use m_helper

use m_constants

implicit none

private :: s_compute_image_points, &
Expand All @@ -39,7 +41,6 @@
type(ghost_point), dimension(:), allocatable :: inner_points
!$acc declare create(ghost_points, inner_points)

integer :: gp_layers !< Number of ghost point layers
integer :: num_gps !< Number of ghost points
integer :: num_inner_gps !< Number of ghost points
!$acc declare create(gp_layers, num_gps, num_inner_gps)
Expand All @@ -49,8 +50,6 @@
!> Allocates memory for the variables in the IBM module
subroutine s_initialize_ibm_module()

gp_layers = 3

if (p > 0) then
@:ALLOCATE(ib_markers%sf(-gp_layers:m+gp_layers, &
-gp_layers:n+gp_layers, -gp_layers:p+gp_layers))
Expand All @@ -69,9 +68,8 @@

@:ACC_SETUP_SFs(ib_markers)
@:ACC_SETUP_SFs(levelset)
! @:ALLOCATE(ib_markers%sf(0:m, 0:n, 0:p))

!$acc enter data copyin(gp_layers, num_gps, num_inner_gps)
!$acc enter data copyin(num_gps, num_inner_gps)

end subroutine s_initialize_ibm_module

Expand Down Expand Up @@ -209,7 +207,9 @@

! Calculate velocity of ghost cell
if (gp%slip) then
norm = levelset_norm%sf(j, k, l, patch_id, :)
norm = levelset_norm%sf(gp%loc(1), gp%loc(2), gp%loc(3), gp%ib_patch_id, :)
buf = sqrt(sum(norm**2))
Copy link
Contributor

@haochey haochey Jan 15, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think these lines are redundant because levelset_norm%sf(gp%loc(1), gp%loc(2), gp%loc(3), gp%ib_patch_id, :) is already a norm vector.

And j = gp%loc(1), k = gp%loc(2) , l = gp%loc(3) , so norm = levelset_norm%sf(j, k, l, patch_id, :) should be fine.

Copy link
Member Author

@sbryngelson sbryngelson Jan 15, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I agree, but what about the other part (line 211 and 222 that have buf = ... and norm = norm/buf? that seems different/new

Copy link
Contributor

@haochey haochey Jan 15, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These two lines seem to normalize the norm to be unit vectors while norm are already unit vectors because levelset_norms' are unit vectors.

Copy link
Member Author

@sbryngelson sbryngelson Jan 16, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i suspect this was not added by @henryleberre for no reason, is there any way to ensure levelset_norm subroutines always normalize without just doing it here?

Copy link
Contributor

@haochey haochey Jan 16, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Make sense to me! It will rely on levelset subroutines to normalize them before output them. I think this is the reason why Henry added these two lines.

Copy link
Contributor

@anandrdbz anandrdbz Jan 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Btw there's a corner case that needs to be taken care of, when norm = 0

Copy link
Member Author

@sbryngelson sbryngelson Jan 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what do you suggest doing in this case, @anandrdbz? You can write code to my branch/PR here.

Copy link
Contributor

@anandrdbz anandrdbz Jan 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay I'll take a look, do I have push permissions ?

Copy link
Member Author

@sbryngelson sbryngelson Jan 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe so.

norm = norm/buf

Check warning on line 212 in src/simulation/m_ibm.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_ibm.fpp#L210-L212 

Added lines #L210 - L212 were not covered by tests
vel_norm_IP = sum(vel_IP*norm)*norm
vel_g = vel_IP - vel_norm_IP
else
Expand Down Expand Up @@ -417,6 +417,7 @@
integer :: i, j, k, l, q !< Iterator variables

num_gps = 0
num_inner_gps = 0

do i = 0, m
do j = 0, n
Expand Down
12 changes: 0 additions & 12 deletions src/simulation/m_weno.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,6 @@ module m_weno
!! of the characteristic decomposition are stored in custom-constructed WENO-
!! stencils (WS) that are annexed to each position of a given scalar field.
!> @{

real(wp), allocatable, dimension(:, :, :, :) :: v_rs_ws_x, v_rs_ws_y, v_rs_ws_z
!> @}

Expand All @@ -52,52 +51,41 @@ module m_weno
!! second dimension identifies the position of its coefficients and the last
!! dimension denotes the cell-location in the relevant coordinate direction.
!> @{

real(wp), target, allocatable, dimension(:, :, :) :: poly_coef_cbL_x
real(wp), target, allocatable, dimension(:, :, :) :: poly_coef_cbL_y
real(wp), target, allocatable, dimension(:, :, :) :: poly_coef_cbL_z

real(wp), target, allocatable, dimension(:, :, :) :: poly_coef_cbR_x
real(wp), target, allocatable, dimension(:, :, :) :: poly_coef_cbR_y
real(wp), target, allocatable, dimension(:, :, :) :: poly_coef_cbR_z

! real(wp), pointer, dimension(:, :, :) :: poly_coef_L => null()
! real(wp), pointer, dimension(:, :, :) :: poly_coef_R => null()
!> @}

!> @name The ideal weights at the left and the right cell-boundaries and at the
!! left and the right quadrature points, in x-, y- and z-directions. Note
!! that the first dimension of the array identifies the weight, while the
!! last denotes the cell-location in the relevant coordinate direction.
!> @{

real(wp), target, allocatable, dimension(:, :) :: d_cbL_x
real(wp), target, allocatable, dimension(:, :) :: d_cbL_y
real(wp), target, allocatable, dimension(:, :) :: d_cbL_z

real(wp), target, allocatable, dimension(:, :) :: d_cbR_x
real(wp), target, allocatable, dimension(:, :) :: d_cbR_y
real(wp), target, allocatable, dimension(:, :) :: d_cbR_z
! real(wp), pointer, dimension(:, :) :: d_L => null()
! real(wp), pointer, dimension(:, :) :: d_R => null()
!> @}

!> @name Smoothness indicator coefficients in the x-, y-, and z-directions. Note
!! that the first array dimension identifies the smoothness indicator, the
!! second identifies the position of its coefficients and the last denotes
!! the cell-location in the relevant coordinate direction.
!> @{

real(wp), target, allocatable, dimension(:, :, :) :: beta_coef_x
real(wp), target, allocatable, dimension(:, :, :) :: beta_coef_y
real(wp), target, allocatable, dimension(:, :, :) :: beta_coef_z
! real(wp), pointer, dimension(:, :, :) :: beta_coef => null()
!> @}

! END: WENO Coefficients

integer :: v_size !< Number of WENO-reconstructed cell-average variables

!$acc declare create(v_size)

!> @name Indical bounds in the s1-, s2- and s3-directions
Expand Down
Loading