Skip to content

Commit

Permalink
Moving Levelset Calculation to Pre-Process (#662)
Browse files Browse the repository at this point in the history
  • Loading branch information
@anshgupta1234
anshgupta1234 authored Nov 12, 2024
1 parent efc9d67 commit daf3acd
Show file tree
Hide file tree
Showing 35 changed files with 1,122 additions and 902 deletions.
20 changes: 20 additions & 0 deletions src/common/m_derived_types.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,16 @@ module m_derived_types
integer, pointer, dimension(:, :, :) :: sf => null()
end type integer_field

!> Derived type for levelset
type levelset_field
real(kind(0d0)), pointer, dimension(:, :, :, :) :: sf => null()
end type levelset_field

!> Derived type for levelset norm
type levelset_norm_field
real(kind(0d0)), pointer, dimension(:, :, :, :, :) :: sf => null()
end type levelset_norm_field

type mpi_io_var
integer, allocatable, dimension(:) :: view
type(scalar_field), allocatable, dimension(:) :: var
Expand All @@ -44,6 +54,16 @@ module m_derived_types
type(integer_field) :: var
end type mpi_io_ib_var

type mpi_io_levelset_var
integer :: view
type(levelset_field) :: var
end type mpi_io_levelset_var

type mpi_io_levelset_norm_var
integer :: view
type(levelset_norm_field) :: var
end type mpi_io_levelset_norm_var

!> Derived type annexing a vector field (VF)
type vector_field
type(scalar_field), allocatable, dimension(:) :: vf !< Vector field
Expand Down
29 changes: 28 additions & 1 deletion src/common/m_mpi_common.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,9 @@ contains
!! @param q_cons_vf Conservative variables
!! @param ib_markers track if a cell is within the immersed boundary
subroutine s_initialize_mpi_data(q_cons_vf, ib_markers)
!! @param levelset closest distance from every cell to the IB
!! @param levelset_norm normalized vector from every cell to the closest point to the IB
subroutine s_initialize_mpi_data(q_cons_vf, ib_markers, levelset, levelset_norm)
type(scalar_field), &
dimension(sys_size), &
Expand All @@ -72,6 +74,14 @@ contains
optional, &
intent(in) :: ib_markers
type(levelset_field), &
optional, &
intent(IN) :: levelset
type(levelset_norm_field), &
optional, &
intent(IN) :: levelset_norm
integer, dimension(num_dims) :: sizes_glb, sizes_loc
integer, dimension(1) :: airfoil_glb, airfoil_loc, airfoil_start
Expand Down Expand Up @@ -137,13 +147,30 @@ contains
#ifdef MFC_PRE_PROCESS
MPI_IO_IB_DATA%var%sf => ib_markers%sf
MPI_IO_levelset_DATA%var%sf => levelset%sf
MPI_IO_levelsetnorm_DATA%var%sf => levelset_norm%sf
#else
MPI_IO_IB_DATA%var%sf => ib_markers%sf(0:m, 0:n, 0:p)
#ifndef MFC_POST_PROCESS
MPI_IO_levelset_DATA%var%sf => levelset%sf(0:m, 0:n, 0:p, 1:num_ibs)
MPI_IO_levelsetnorm_DATA%var%sf => levelset_norm%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3)
#endif
#endif
call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, &
MPI_ORDER_FORTRAN, MPI_INTEGER, MPI_IO_IB_DATA%view, ierr)
call MPI_TYPE_COMMIT(MPI_IO_IB_DATA%view, ierr)
#ifndef MFC_POST_PROCESS
call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, &
MPI_ORDER_FORTRAN, MPI_DOUBLE_PRECISION, MPI_IO_levelset_DATA%view, ierr)
call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, &
MPI_ORDER_FORTRAN, MPI_DOUBLE_PRECISION, MPI_IO_levelsetnorm_DATA%view, ierr)
call MPI_TYPE_COMMIT(MPI_IO_levelset_DATA%view, ierr)
call MPI_TYPE_COMMIT(MPI_IO_levelsetnorm_DATA%view, ierr)
#endif
end if
#ifndef MFC_POST_PROCESS
Expand Down
7 changes: 7 additions & 0 deletions src/post_process/m_global_parameters.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,10 +155,14 @@ module m_global_parameters
integer, allocatable, dimension(:) :: start_idx !<
!! Starting cell-center index of local processor in global grid

integer :: num_ibs !< Number of immersed boundaries

#ifdef MFC_MPI

type(mpi_io_var), public :: MPI_IO_DATA
type(mpi_io_ib_var), public :: MPI_IO_IB_DATA
type(mpi_io_levelset_var), public :: MPI_IO_levelset_DATA
type(mpi_io_levelset_norm_var), public :: MPI_IO_levelsetnorm_DATA

#endif

Expand Down Expand Up @@ -401,6 +405,9 @@ contains
surface_tension = .false.
adv_n = .false.

! IBM
num_ibs = dflt_int

end subroutine s_assign_default_values_to_user_inputs

!> Computation of parameters, allocation procedures, and/or
Expand Down
2 changes: 1 addition & 1 deletion src/post_process/m_mpi_proxy.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,7 @@ contains
& 'model_eqns', 'num_fluids', 'bc_x%beg', 'bc_x%end', 'bc_y%beg', &
& 'bc_y%end', 'bc_z%beg', 'bc_z%end', 'flux_lim', 'format', &
& 'precision', 'fd_order', 'thermal', 'nb', 'relax_model', &
& 'n_start' ]
& 'n_start', 'num_ibs' ]
call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
#:endfor
Expand Down
2 changes: 1 addition & 1 deletion src/post_process/m_start_up.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ subroutine s_read_input_file
parallel_io, rhoref, pref, bubbles, qbmm, sigR, &
R0ref, nb, polytropic, thermal, Ca, Web, Re_inv, &
polydisperse, poly_sigma, file_per_process, relax, &
relax_model, cf_wrt, sigma, adv_n, ib, &
relax_model, cf_wrt, sigma, adv_n, ib, num_ibs, &
cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, &
cfl_target, surface_tension

Expand Down
111 changes: 55 additions & 56 deletions src/simulation/m_compute_levelset.fpp → src/pre_process/m_compute_levelset.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,12 +38,11 @@ module m_compute_levelset

contains

!> Initialize IBM module
subroutine s_compute_circle_levelset(levelset, levelset_norm, ib_patch_id)

real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
integer, intent(in) :: ib_patch_id
type(levelset_field), intent(INOUT) :: levelset
type(levelset_norm_field), intent(INOUT) :: levelset_norm
integer, intent(IN) :: ib_patch_id

real(kind(0d0)) :: radius, dist
real(kind(0d0)) :: x_centroid, y_centroid
Expand All @@ -62,11 +61,11 @@ contains
dist_vec(2) = y_cc(j) - y_centroid
dist_vec(3) = 0
dist = dsqrt(sum(dist_vec**2))
levelset(i, j, 0, ib_patch_id) = dist - radius
levelset%sf(i, j, 0, ib_patch_id) = dist - radius
if (dist == 0) then
levelset_norm(i, j, 0, ib_patch_id, :) = 0
levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0
else
levelset_norm(i, j, 0, ib_patch_id, :) = &
levelset_norm%sf(i, j, 0, ib_patch_id, :) = &
dist_vec(:)/dist
end if

Expand All @@ -77,9 +76,9 @@ contains

subroutine s_compute_airfoil_levelset(levelset, levelset_norm, ib_patch_id)

real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
integer, intent(in) :: ib_patch_id
type(levelset_field), intent(INOUT) :: levelset
type(levelset_norm_field), intent(INOUT) :: levelset_norm
integer, intent(IN) :: ib_patch_id

real(kind(0d0)) :: radius, dist, global_dist
integer :: global_id
Expand Down Expand Up @@ -145,11 +144,11 @@ contains
dist = global_dist
end if

levelset(i, j, 0, ib_patch_id) = dist
levelset%sf(i, j, 0, ib_patch_id) = dist
if (dist == 0) then
levelset_norm(i, j, 0, ib_patch_id, :) = 0
levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0
else
levelset_norm(i, j, 0, ib_patch_id, :) = &
levelset_norm%sf(i, j, 0, ib_patch_id, :) = &
dist_vec(:)/dist
end if

Expand All @@ -160,9 +159,9 @@ contains

subroutine s_compute_3D_airfoil_levelset(levelset, levelset_norm, ib_patch_id)

real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
integer, intent(in) :: ib_patch_id
type(levelset_field), intent(INOUT) :: levelset
type(levelset_norm_field), intent(INOUT) :: levelset_norm
integer, intent(IN) :: ib_patch_id

real(kind(0d0)) :: radius, dist, dist_surf, dist_side, global_dist
integer :: global_id
Expand Down Expand Up @@ -237,18 +236,18 @@ contains
dist_side = min(abs(z_cc(l) - z_min), abs(z_max - z_cc(l)))

if (dist_side < dist_surf) then
levelset(i, j, l, ib_patch_id) = dist_side
levelset%sf(i, j, l, ib_patch_id) = dist_side
if (dist_side == abs(z_cc(l) - z_min)) then
levelset_norm(i, j, l, ib_patch_id, :) = (/0, 0, -1/)
levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, -1/)
else
levelset_norm(i, j, l, ib_patch_id, :) = (/0, 0, 1/)
levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, 1/)
end if
else
levelset(i, j, l, ib_patch_id) = dist_surf
levelset%sf(i, j, l, ib_patch_id) = dist_surf
if (dist == 0) then
levelset_norm(i, j, l, ib_patch_id, :) = 0
levelset_norm%sf(i, j, l, ib_patch_id, :) = 0
else
levelset_norm(i, j, l, ib_patch_id, :) = &
levelset_norm%sf(i, j, l, ib_patch_id, :) = &
dist_vec(:)/dist_surf
end if
end if
Expand All @@ -262,10 +261,10 @@ contains
!> Initialize IBM module
subroutine s_compute_rectangle_levelset(levelset, levelset_norm, ib_patch_id)

real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
integer, intent(in) :: ib_patch_id
type(levelset_field), intent(INOUT) :: levelset
type(levelset_norm_field), intent(INOUT) :: levelset_norm

integer :: ib_patch_id
real(kind(0d0)) :: top_right(2), bottom_left(2)
real(kind(0d0)) :: x, y, min_dist
real(kind(0d0)) :: side_dists(4)
Expand Down Expand Up @@ -303,37 +302,37 @@ contains
end if

if (min_dist == abs(side_dists(1))) then
levelset(i, j, 0, ib_patch_id) = side_dists(1)
levelset%sf(i, j, 0, ib_patch_id) = side_dists(1)
if (side_dists(1) == 0) then
levelset_norm(i, j, 0, ib_patch_id, 1) = 0d0
levelset_norm%sf(i, j, 0, ib_patch_id, 1) = 0d0
else
levelset_norm(i, j, 0, ib_patch_id, 1) = side_dists(1)/ &
abs(side_dists(1))
levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(1)/ &
abs(side_dists(1))
end if

else if (min_dist == abs(side_dists(2))) then
levelset(i, j, 0, ib_patch_id) = side_dists(2)
levelset%sf(i, j, 0, ib_patch_id) = side_dists(2)
if (side_dists(2) == 0) then
levelset_norm(i, j, 0, ib_patch_id, 1) = 0d0
levelset_norm%sf(i, j, 0, ib_patch_id, 1) = 0d0
else
levelset_norm(i, j, 0, ib_patch_id, 1) = side_dists(2)/ &
abs(side_dists(2))
levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(2)/ &
abs(side_dists(2))
end if

else if (min_dist == abs(side_dists(3))) then
if (side_dists(3) == 0) then
levelset_norm(i, j, 0, ib_patch_id, 1) = 0d0
levelset_norm%sf(i, j, 0, ib_patch_id, 1) = 0d0
else
levelset_norm(i, j, 0, ib_patch_id, 1) = side_dists(3)/ &
abs(side_dists(3))
levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(3)/ &
abs(side_dists(3))
end if

else if (min_dist == abs(side_dists(4))) then
if (side_dists(4) == 0) then
levelset_norm(i, j, 0, ib_patch_id, 1) = 0d0
levelset_norm%sf(i, j, 0, ib_patch_id, 1) = 0d0
else
levelset_norm(i, j, 0, ib_patch_id, 1) = side_dists(4)/ &
abs(side_dists(4))
levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(4)/ &
abs(side_dists(4))
end if

end if
Expand All @@ -347,9 +346,9 @@ contains

subroutine s_compute_sphere_levelset(levelset, levelset_norm, ib_patch_id)

real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
integer, intent(in) :: ib_patch_id
type(levelset_field), intent(INOUT) :: levelset
type(levelset_norm_field), intent(INOUT) :: levelset_norm
integer, intent(IN) :: ib_patch_id

real(kind(0d0)) :: radius, dist
real(kind(0d0)) :: x_centroid, y_centroid, z_centroid
Expand All @@ -369,11 +368,11 @@ contains
dist_vec(2) = y_cc(j) - y_centroid
dist_vec(3) = z_cc(k) - z_centroid
dist = dsqrt(sum(dist_vec**2))
levelset(i, j, k, ib_patch_id) = dist - radius
levelset%sf(i, j, k, ib_patch_id) = dist - radius
if (dist == 0) then
levelset_norm(i, j, k, ib_patch_id, :) = (/1, 0, 0/)
levelset_norm%sf(i, j, k, ib_patch_id, :) = (/1, 0, 0/)
else
levelset_norm(i, j, k, ib_patch_id, :) = &
levelset_norm%sf(i, j, k, ib_patch_id, :) = &
dist_vec(:)/dist
end if
end do
Expand All @@ -384,9 +383,9 @@ contains

subroutine s_compute_cylinder_levelset(levelset, levelset_norm, ib_patch_id)

real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
integer, intent(in) :: ib_patch_id
type(levelset_field), intent(INOUT) :: levelset
type(levelset_norm_field), intent(INOUT) :: levelset_norm
integer, intent(IN) :: ib_patch_id

real(kind(0d0)) :: radius, dist
real(kind(0d0)) :: x_centroid, y_centroid, z_centroid
Expand Down Expand Up @@ -434,20 +433,20 @@ contains
- radius

if (dist_side < abs(dist_surface)) then
levelset(i, j, k, ib_patch_id) = -dist_side
levelset%sf(i, j, k, ib_patch_id) = -dist_side
if (dist_side == abs(side_pos - boundary%beg)) then
levelset_norm(i, j, k, ib_patch_id, :) = -dist_sides_vec
levelset_norm%sf(i, j, k, ib_patch_id, :) = -dist_sides_vec
else
levelset_norm(i, j, k, ib_patch_id, :) = dist_sides_vec
levelset_norm%sf(i, j, k, ib_patch_id, :) = dist_sides_vec
end if
else
levelset(i, j, k, ib_patch_id) = dist_surface
levelset%sf(i, j, k, ib_patch_id) = dist_surface

levelset_norm(i, j, k, ib_patch_id, :) = &
levelset_norm%sf(i, j, k, ib_patch_id, :) = &
(pos_vec - centroid_vec)*dist_surface_vec
levelset_norm(i, j, k, ib_patch_id, :) = &
levelset_norm(i, j, k, ib_patch_id, :)/ &
norm2(levelset_norm(i, j, k, ib_patch_id, :))
levelset_norm%sf(i, j, k, ib_patch_id, :) = &
levelset_norm%sf(i, j, k, ib_patch_id, :)/ &
norm2(levelset_norm%sf(i, j, k, ib_patch_id, :))
end if
end do
end do
Expand Down
Loading

0 comments on commit daf3acd

Please sign in to comment.