diff --git a/.vscode/settings.json b/.vscode/settings.json index ab9cd6f770..5c880f8ded 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -49,7 +49,7 @@ // Enable ONLY fortls language server "fortran.enableLanguageServer": true, "fortran.languageServer": "fortls", - "fortran.fortls.disabled": false, + "fortran.fortls.disabled": true, "fortran.fortls.path": "fortls", // Try to disable any built-in language features diff --git a/CMakeLists.txt b/CMakeLists.txt index 0b349eb394..1b6970ad9f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -328,6 +328,13 @@ macro(HANDLE_SOURCES target useCommon) "${CMAKE_BINARY_DIR}/modules/${target}/*.fpp") if (${useCommon}) file(GLOB common_FPPs CONFIGURE_DEPENDS "${common_DIR}/*.fpp") + + # If we're building post_process, exclude m_compute_levelset.fpp + if("${target}" STREQUAL "post_process") + list(FILTER common_FPPs EXCLUDE REGEX ".*/m_compute_levelset\.fpp$") + list(FILTER common_FPPs EXCLUDE REGEX ".*/m_ib_patches\.fpp$") + endif() + list(APPEND ${target}_FPPs ${common_FPPs}) endif() @@ -443,7 +450,7 @@ function(MFC_SETUP_TARGET) if (ARGS_SILO) find_package(SILO REQUIRED) - target_link_libraries(${a_target} PRIVATE SILO::SILO) + target_link_libraries(${a_target} PRIVATE SILO::SILO stdc++) endif() if (ARGS_HDF5) diff --git a/src/pre_process/include/1dHardcodedIC.fpp b/src/common/include/1dHardcodedIC.fpp similarity index 100% rename from src/pre_process/include/1dHardcodedIC.fpp rename to src/common/include/1dHardcodedIC.fpp diff --git a/src/pre_process/include/2dHardcodedIC.fpp b/src/common/include/2dHardcodedIC.fpp similarity index 100% rename from src/pre_process/include/2dHardcodedIC.fpp rename to src/common/include/2dHardcodedIC.fpp diff --git a/src/pre_process/include/3dHardcodedIC.fpp b/src/common/include/3dHardcodedIC.fpp similarity index 100% rename from src/pre_process/include/3dHardcodedIC.fpp rename to src/common/include/3dHardcodedIC.fpp diff --git a/src/pre_process/include/ExtrusionHardcodedIC.fpp b/src/common/include/ExtrusionHardcodedIC.fpp similarity index 100% rename from src/pre_process/include/ExtrusionHardcodedIC.fpp rename to src/common/include/ExtrusionHardcodedIC.fpp diff --git a/src/common/include/case.fpp b/src/common/include/case.fpp index ad2e0b1a94..84d1d2cc14 100644 --- a/src/common/include/case.fpp +++ b/src/common/include/case.fpp @@ -4,4 +4,5 @@ ! For pre-process. #:def analytical() + #:enddef diff --git a/src/pre_process/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp similarity index 100% rename from src/pre_process/m_compute_levelset.fpp rename to src/common/m_compute_levelset.fpp diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 38674af615..9c55e8b895 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -318,6 +318,11 @@ module m_derived_types real(wp) :: model_threshold !< !! Threshold to turn on smoothen STL patch. + + !! Patch conditions for moving imersed boundaries + integer :: moving_ibm ! 0 for no moving, 1 for moving, 2 for moving on forced path + real(wp), dimension(1:3) :: vel + end type ib_patch_parameters !> Derived type annexing the physical parameters (PP) of the fluids. These diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp new file mode 100644 index 0000000000..3403b4c7e5 --- /dev/null +++ b/src/common/m_ib_patches.fpp @@ -0,0 +1,1040 @@ +!> +!! @file m_patches.fpp +!! @brief Contains module m_patches + +#:include 'case.fpp' +#:include 'ExtrusionHardcodedIC.fpp' +#:include '1dHardcodedIC.fpp' +#:include '2dHardcodedIC.fpp' +#:include '3dHardcodedIC.fpp' +#:include 'macros.fpp' + +module m_ib_patches + + use m_model ! Subroutine(s) related to STL files + + use m_derived_types ! Definitions of the derived types + + use m_global_parameters !< Definitions of the global parameters + + use m_helper_basic !< Functions to compare floating point numbers + + use m_helper + + use m_compute_levelset ! Subroutines to calculate levelsets for IBs + + use m_mpi_common + + implicit none + + private; public :: s_apply_ib_patches + + real(wp) :: x_centroid, y_centroid, z_centroid + real(wp) :: length_x, length_y, length_z + + integer :: smooth_patch_id + real(wp) :: smooth_coeff !< + !! These variables are analogous in both meaning and use to the similarly + !! named components in the ic_patch_parameters type (see m_derived_types.f90 + !! for additional details). They are employed as a means to more concisely + !! perform the actions necessary to lay out a particular patch on the grid. + + real(wp) :: eta !< + !! In the case that smoothing of patch boundaries is enabled and the boundary + !! between two adjacent patches is to be smeared out, this variable's purpose + !! is to act as a pseudo volume fraction to indicate the contribution of each + !! patch toward the composition of a cell's fluid state. + + real(wp) :: cart_x, cart_y, cart_z + real(wp) :: sph_phi !< + !! Variables to be used to hold cell locations in Cartesian coordinates if + !! 3D simulation is using cylindrical coordinates + + type(bounds_info) :: x_boundary, y_boundary, z_boundary !< + !! These variables combine the centroid and length parameters associated with + !! a particular patch to yield the locations of the patch boundaries in the + !! x-, y- and z-coordinate directions. They are used as a means to concisely + !! perform the actions necessary to lay out a particular patch on the grid. + + character(len=5) :: istr ! string to store int to string result for error checking + +contains + + impure subroutine s_apply_ib_patches(ib_markers_sf, levelset, levelset_norm) + + integer, dimension(:, :, :), intent(inout), optional :: ib_markers_sf + type(levelset_field), intent(inout), optional :: levelset !< Levelset determined by models + type(levelset_norm_field), intent(inout), optional :: levelset_norm !< Levelset_norm determined by models + + integer :: i + + ! 3D Patch Geometries + if (p > 0) then + + !> IB Patches + !> @{ + ! Spherical patch + do i = 1, num_ibs + + if (patch_ib(i)%geometry == 8) then + call s_ib_sphere(i, ib_markers_sf) + call s_sphere_levelset(i, levelset, levelset_norm) + elseif (patch_ib(i)%geometry == 9) then + call s_ib_cuboid(i, ib_markers_sf) + call s_cuboid_levelset(i, levelset, levelset_norm) + elseif (patch_ib(i)%geometry == 10) then + call s_ib_cylinder(i, ib_markers_sf) + call s_cylinder_levelset(i, levelset, levelset_norm) + elseif (patch_ib(i)%geometry == 11) then +#ifdef MFC_PRE_PROCESS + call s_ib_3D_airfoil(i, ib_markers_sf) + call s_3D_airfoil_levelset(i, levelset, levelset_norm) +#endif + ! STL+IBM patch + elseif (patch_ib(i)%geometry == 12) then + call s_ib_model(i, ib_markers_sf, levelset, levelset_norm) + end if + end do + !> @} + + ! 2D Patch Geometries + elseif (n > 0) then + + !> IB Patches + !> @{ + do i = 1, num_ibs + if (patch_ib(i)%geometry == 2) then + call s_ib_circle(i, ib_markers_sf) + call s_circle_levelset(i, levelset, levelset_norm) + elseif (patch_ib(i)%geometry == 3) then + call s_ib_rectangle(i, ib_markers_sf) + call s_rectangle_levelset(i, levelset, levelset_norm) + elseif (patch_ib(i)%geometry == 4) then +#ifdef MFC_PRE_PROCESS + call s_ib_airfoil(i, ib_markers_sf) + call s_airfoil_levelset(i, levelset, levelset_norm) +#endif + ! STL+IBM patch + elseif (patch_ib(i)%geometry == 5) then + call s_ib_model(i, ib_markers_sf, levelset, levelset_norm) + end if + end do + !> @} + + end if + + end subroutine s_apply_ib_patches + + !> The circular patch is a 2D geometry that may be used, for + !! example, in creating a bubble or a droplet. The geometry + !! of the patch is well-defined when its centroid and radius + !! are provided. Note that the circular patch DOES allow for + !! the smoothing of its boundary. + !! @param patch_id is the patch identifier + !! @param ib_markers_sf Array to track patch ids + !! @param ib True if this patch is an immersed boundary + subroutine s_ib_circle(patch_id, ib_markers_sf) + + integer, intent(in) :: patch_id + integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + + real(wp) :: radius + + integer :: i, j, k !< Generic loop iterators + + ! Transferring the circular patch's radius, centroid, smearing patch + ! identity and smearing coefficient information + + x_centroid = patch_ib(patch_id)%x_centroid + y_centroid = patch_ib(patch_id)%y_centroid + radius = patch_ib(patch_id)%radius + + ! Initializing the pseudo volume fraction value to 1. The value will + ! be modified as the patch is laid out on the grid, but only in the + ! case that smoothing of the circular patch's boundary is enabled. + eta = 1._wp + + ! Checking whether the circle covers a particular cell in the domain + ! and verifying whether the current patch has permission to write to + ! that cell. If both queries check out, the primitive variables of + ! the current patch are assigned to this cell. + + ! TODO :: THIS SETS ib_markers_sf TO HOLD THE PATCH ID, BUT WE NEED TO + ! NOW ALSO SEARCH FOR OTHER POINTS TO DELETE THE CURRENT PATCH ID + + do j = 0, n + do i = 0, m + if ((x_cc(i) - x_centroid)**2 & + + (y_cc(j) - y_centroid)**2 <= radius**2) & + then + ib_markers_sf(i, j, 0) = patch_id + end if + end do + end do + + end subroutine s_ib_circle + + !! @param patch_id is the patch identifier + !! @param ib_markers_sf Array to track patch ids + subroutine s_ib_airfoil(patch_id, ib_markers_sf) + + integer, intent(in) :: patch_id + integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + + real(wp) :: x0, y0, f, x_act, y_act, ca_in, pa, ma, ta, theta + real(wp) :: xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c + integer :: i, j, k + integer :: Np1, Np2 + + x0 = patch_ib(patch_id)%x_centroid + y0 = patch_ib(patch_id)%y_centroid + ca_in = patch_ib(patch_id)%c + pa = patch_ib(patch_id)%p + ma = patch_ib(patch_id)%m + ta = patch_ib(patch_id)%t + theta = pi*patch_ib(patch_id)%theta/180._wp + + ! rank(dx) is not consistent between pre_process and simulation. This IFDEF prevents compilation errors +#ifdef MFC_PRE_PROCESS + Np1 = int((pa*ca_in/dx)*20) + Np2 = int(((ca_in - pa*ca_in)/dx)*20) +#else + Np1 = int((pa*ca_in/dx(0))*20) + Np2 = int(((ca_in - pa*ca_in)/dx(0))*20) +#endif + Np = Np1 + Np2 + 1 + + allocate (airfoil_grid_u(1:Np)) + allocate (airfoil_grid_l(1:Np)) + + airfoil_grid_u(1)%x = x0 + airfoil_grid_u(1)%y = y0 + + airfoil_grid_l(1)%x = x0 + airfoil_grid_l(1)%y = y0 + + eta = 1._wp + + do i = 1, Np1 + Np2 - 1 + if (i <= Np1) then + xc = x0 + i*(pa*ca_in/Np1) + xa = (xc - x0)/ca_in + yc = (ma/pa**2)*(2*pa*xa - xa**2) + dycdxc = (2*ma/pa**2)*(pa - xa) + else + xc = x0 + pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2) + xa = (xc - x0)/ca_in + yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) + dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) + end if + + yt = (5._wp*ta)*(0.2969_wp*xa**0.5_wp - 0.126_wp*xa - 0.3516_wp*xa**2._wp + 0.2843_wp*xa**3 - 0.1015_wp*xa**4) + sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp + cos_c = 1/(1 + dycdxc**2)**0.5_wp + + xu = xa - yt*sin_c + yu = yc + yt*cos_c + + xl = xa + yt*sin_c + yl = yc - yt*cos_c + + xu = xu*ca_in + x0 + yu = yu*ca_in + y0 + + xl = xl*ca_in + x0 + yl = yl*ca_in + y0 + + airfoil_grid_u(i + 1)%x = xu + airfoil_grid_u(i + 1)%y = yu + + airfoil_grid_l(i + 1)%x = xl + airfoil_grid_l(i + 1)%y = yl + + end do + + airfoil_grid_u(Np)%x = x0 + ca_in + airfoil_grid_u(Np)%y = y0 + + airfoil_grid_l(Np)%x = x0 + ca_in + airfoil_grid_l(Np)%y = y0 + + do j = 0, n + do i = 0, m + + if (.not. f_is_default(patch_ib(patch_id)%theta)) then + x_act = (x_cc(i) - x0)*cos(theta) - (y_cc(j) - y0)*sin(theta) + x0 + y_act = (x_cc(i) - x0)*sin(theta) + (y_cc(j) - y0)*cos(theta) + y0 + else + x_act = x_cc(i) + y_act = y_cc(j) + end if + + if (x_act >= x0 .and. x_act <= x0 + ca_in) then + xa = (x_act - x0)/ca_in + if (xa <= pa) then + yc = (ma/pa**2)*(2*pa*xa - xa**2) + dycdxc = (2*ma/pa**2)*(pa - xa) + else + yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) + dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) + end if + if (y_act >= y0) then + k = 1 + do while (airfoil_grid_u(k)%x < x_act) + k = k + 1 + end do + if (f_approx_equal(airfoil_grid_u(k)%x, x_act)) then + if (y_act <= airfoil_grid_u(k)%y) then + !!IB + ib_markers_sf(i, j, 0) = patch_id + end if + else + f = (airfoil_grid_u(k)%x - x_act)/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x) + if (y_act <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then + !!IB + ib_markers_sf(i, j, 0) = patch_id + end if + end if + else + k = 1 + do while (airfoil_grid_l(k)%x < x_act) + k = k + 1 + end do + if (f_approx_equal(airfoil_grid_l(k)%x, x_act)) then + if (y_act >= airfoil_grid_l(k)%y) then + !!IB + ib_markers_sf(i, j, 0) = patch_id + end if + else + f = (airfoil_grid_l(k)%x - x_act)/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x) + + if (y_act >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then + !!IB + ib_markers_sf(i, j, 0) = patch_id + end if + end if + end if + end if + end do + end do + + if (.not. f_is_default(patch_ib(patch_id)%theta)) then + do i = 1, Np + airfoil_grid_l(i)%x = (airfoil_grid_l(i)%x - x0)*cos(theta) + (airfoil_grid_l(i)%y - y0)*sin(theta) + x0 + airfoil_grid_l(i)%y = -1._wp*(airfoil_grid_l(i)%x - x0)*sin(theta) + (airfoil_grid_l(i)%y - y0)*cos(theta) + y0 + + airfoil_grid_u(i)%x = (airfoil_grid_u(i)%x - x0)*cos(theta) + (airfoil_grid_u(i)%y - y0)*sin(theta) + x0 + airfoil_grid_u(i)%y = -1._wp*(airfoil_grid_u(i)%x - x0)*sin(theta) + (airfoil_grid_u(i)%y - y0)*cos(theta) + y0 + end do + end if + + end subroutine s_ib_airfoil + + !! @param patch_id is the patch identifier + !! @param ib_markers_sf Array to track patch ids + !! @param ib True if this patch is an immersed boundary + subroutine s_ib_3D_airfoil(patch_id, ib_markers_sf) + + integer, intent(in) :: patch_id + integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + + real(wp) :: x0, y0, z0, lz, z_max, z_min, f, x_act, y_act, ca_in, pa, ma, ta, theta, xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c + integer :: i, j, k, l + integer :: Np1, Np2 + + x0 = patch_ib(patch_id)%x_centroid + y0 = patch_ib(patch_id)%y_centroid + z0 = patch_ib(patch_id)%z_centroid + lz = patch_ib(patch_id)%length_z + ca_in = patch_ib(patch_id)%c + pa = patch_ib(patch_id)%p + ma = patch_ib(patch_id)%m + ta = patch_ib(patch_id)%t + theta = pi*patch_ib(patch_id)%theta/180._wp + + ! rank(dx) is not consistent between pre_process and simulation. This IFDEF prevents compilation errors +#ifdef MFC_PRE_PROCESS + Np1 = int((pa*ca_in/dx)*20) + Np2 = int(((ca_in - pa*ca_in)/dx)*20) +#else + Np1 = int((pa*ca_in/dx(0))*20) + Np2 = int(((ca_in - pa*ca_in)/dx(0))*20) +#endif + Np = Np1 + Np2 + 1 + + allocate (airfoil_grid_u(1:Np)) + allocate (airfoil_grid_l(1:Np)) + + airfoil_grid_u(1)%x = x0 + airfoil_grid_u(1)%y = y0 + + airfoil_grid_l(1)%x = x0 + airfoil_grid_l(1)%y = y0 + + z_max = z0 + lz/2 + z_min = z0 - lz/2 + + eta = 1._wp + + do i = 1, Np1 + Np2 - 1 + if (i <= Np1) then + xc = x0 + i*(pa*ca_in/Np1) + xa = (xc - x0)/ca_in + yc = (ma/pa**2)*(2*pa*xa - xa**2) + dycdxc = (2*ma/pa**2)*(pa - xa) + else + xc = x0 + pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2) + xa = (xc - x0)/ca_in + yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) + dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) + end if + + yt = (5._wp*ta)*(0.2969_wp*xa**0.5_wp - 0.126_wp*xa - 0.3516_wp*xa**2._wp + 0.2843_wp*xa**3 - 0.1015_wp*xa**4) + sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp + cos_c = 1/(1 + dycdxc**2)**0.5_wp + + xu = xa - yt*sin_c + yu = yc + yt*cos_c + + xl = xa + yt*sin_c + yl = yc - yt*cos_c + + xu = xu*ca_in + x0 + yu = yu*ca_in + y0 + + xl = xl*ca_in + x0 + yl = yl*ca_in + y0 + + airfoil_grid_u(i + 1)%x = xu + airfoil_grid_u(i + 1)%y = yu + + airfoil_grid_l(i + 1)%x = xl + airfoil_grid_l(i + 1)%y = yl + + end do + + airfoil_grid_u(Np)%x = x0 + ca_in + airfoil_grid_u(Np)%y = y0 + + airfoil_grid_l(Np)%x = x0 + ca_in + airfoil_grid_l(Np)%y = y0 + + do l = 0, p + if (z_cc(l) >= z_min .and. z_cc(l) <= z_max) then + do j = 0, n + do i = 0, m + + if (.not. f_is_default(patch_ib(patch_id)%theta)) then + x_act = (x_cc(i) - x0)*cos(theta) - (y_cc(j) - y0)*sin(theta) + x0 + y_act = (x_cc(i) - x0)*sin(theta) + (y_cc(j) - y0)*cos(theta) + y0 + else + x_act = x_cc(i) + y_act = y_cc(j) + end if + + if (x_act >= x0 .and. x_act <= x0 + ca_in) then + xa = (x_act - x0)/ca_in + if (xa <= pa) then + yc = (ma/pa**2)*(2*pa*xa - xa**2) + dycdxc = (2*ma/pa**2)*(pa - xa) + else + yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) + dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) + end if + if (y_act >= y0) then + k = 1 + do while (airfoil_grid_u(k)%x < x_act) + k = k + 1 + end do + if (f_approx_equal(airfoil_grid_u(k)%x, x_act)) then + if (y_act <= airfoil_grid_u(k)%y) then + !!IB + ib_markers_sf(i, j, l) = patch_id + end if + else + f = (airfoil_grid_u(k)%x - x_act)/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x) + if (y_act <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then + !!IB + ib_markers_sf(i, j, l) = patch_id + end if + end if + else + k = 1 + do while (airfoil_grid_l(k)%x < x_act) + k = k + 1 + end do + if (f_approx_equal(airfoil_grid_l(k)%x, x_act)) then + if (y_act >= airfoil_grid_l(k)%y) then + !!IB + ib_markers_sf(i, j, l) = patch_id + end if + else + f = (airfoil_grid_l(k)%x - x_act)/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x) + + if (y_act >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then + !!IB + ib_markers_sf(i, j, l) = patch_id + end if + end if + end if + end if + end do + end do + end if + end do + + if (.not. f_is_default(patch_ib(patch_id)%theta)) then + do i = 1, Np + airfoil_grid_l(i)%x = (airfoil_grid_l(i)%x - x0)*cos(theta) + (airfoil_grid_l(i)%y - y0)*sin(theta) + x0 + airfoil_grid_l(i)%y = -1._wp*(airfoil_grid_l(i)%x - x0)*sin(theta) + (airfoil_grid_l(i)%y - y0)*cos(theta) + y0 + + airfoil_grid_u(i)%x = (airfoil_grid_u(i)%x - x0)*cos(theta) + (airfoil_grid_u(i)%y - y0)*sin(theta) + x0 + airfoil_grid_u(i)%y = -1._wp*(airfoil_grid_u(i)%x - x0)*sin(theta) + (airfoil_grid_u(i)%y - y0)*cos(theta) + y0 + end do + end if + + end subroutine s_ib_3D_airfoil + + !> The rectangular patch is a 2D geometry that may be used, + !! for example, in creating a solid boundary, or pre-/post- + !! shock region, in alignment with the axes of the Cartesian + !! coordinate system. The geometry of such a patch is well- + !! defined when its centroid and lengths in the x- and y- + !! coordinate directions are provided. Please note that the + !! rectangular patch DOES NOT allow for the smoothing of its + !! boundaries. + !! @param patch_id is the patch identifier + !! @param ib_markers_sf Array to track patch ids + !! @param ib True if this patch is an immersed boundary + subroutine s_ib_rectangle(patch_id, ib_markers_sf) + + integer, intent(in) :: patch_id + integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + + integer :: i, j, k !< generic loop iterators + real(wp) :: pi_inf, gamma, lit_gamma !< Equation of state parameters + + pi_inf = fluid_pp(1)%pi_inf + gamma = fluid_pp(1)%gamma + lit_gamma = (1._wp + gamma)/gamma + + ! Transferring the rectangle's centroid and length information + x_centroid = patch_ib(patch_id)%x_centroid + y_centroid = patch_ib(patch_id)%y_centroid + length_x = patch_ib(patch_id)%length_x + length_y = patch_ib(patch_id)%length_y + + ! Computing the beginning and the end x- and y-coordinates of the + ! rectangle based on its centroid and lengths + x_boundary%beg = x_centroid - 0.5_wp*length_x + x_boundary%end = x_centroid + 0.5_wp*length_x + y_boundary%beg = y_centroid - 0.5_wp*length_y + y_boundary%end = y_centroid + 0.5_wp*length_y + + ! Since the rectangular patch does not allow for its boundaries to + ! be smoothed out, the pseudo volume fraction is set to 1 to ensure + ! that only the current patch contributes to the fluid state in the + ! cells that this patch covers. + eta = 1._wp + + ! Checking whether the rectangle covers a particular cell in the + ! domain and verifying whether the current patch has the permission + ! to write to that cell. If both queries check out, the primitive + ! variables of the current patch are assigned to this cell. + do j = 0, n + do i = 0, m + if (x_boundary%beg <= x_cc(i) .and. & + x_boundary%end >= x_cc(i) .and. & + y_boundary%beg <= y_cc(j) .and. & + y_boundary%end >= y_cc(j)) then + + ! Updating the patch identities bookkeeping variable + ib_markers_sf(i, j, 0) = patch_id + + end if + end do + end do + + end subroutine s_ib_rectangle + + !> The spherical patch is a 3D geometry that may be used, + !! for example, in creating a bubble or a droplet. The patch + !! geometry is well-defined when its centroid and radius are + !! provided. Please note that the spherical patch DOES allow + !! for the smoothing of its boundary. + !! @param patch_id is the patch identifier + !! @param ib_markers_sf Array to track patch ids + !! @param ib True if this patch is an immersed boundary + subroutine s_ib_sphere(patch_id, ib_markers_sf) + + integer, intent(in) :: patch_id + integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + + ! Generic loop iterators + integer :: i, j, k + real(wp) :: radius + + !! Variables to initialize the pressure field that corresponds to the + !! bubble-collapse test case found in Tiwari et al. (2013) + + ! Transferring spherical patch's radius, centroid, smoothing patch + ! identity and smoothing coefficient information + x_centroid = patch_ib(patch_id)%x_centroid + y_centroid = patch_ib(patch_id)%y_centroid + z_centroid = patch_ib(patch_id)%z_centroid + radius = patch_ib(patch_id)%radius + + ! Initializing the pseudo volume fraction value to 1. The value will + ! be modified as the patch is laid out on the grid, but only in the + ! case that smoothing of the spherical patch's boundary is enabled. + eta = 1._wp + + ! Checking whether the sphere covers a particular cell in the domain + ! and verifying whether the current patch has permission to write to + ! that cell. If both queries check out, the primitive variables of + ! the current patch are assigned to this cell. + do k = 0, p + do j = 0, n + do i = 0, m + if (grid_geometry == 3) then + call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) + else + cart_y = y_cc(j) + cart_z = z_cc(k) + end if + ! Updating the patch identities bookkeeping variable + if (((x_cc(i) - x_centroid)**2 & + + (cart_y - y_centroid)**2 & + + (cart_z - z_centroid)**2 <= radius**2)) then + ib_markers_sf(i, j, k) = patch_id + end if + end do + end do + end do + + end subroutine s_ib_sphere + + !> The cuboidal patch is a 3D geometry that may be used, for + !! example, in creating a solid boundary, or pre-/post-shock + !! region, which is aligned with the axes of the Cartesian + !! coordinate system. The geometry of such a patch is well- + !! defined when its centroid and lengths in the x-, y- and + !! z-coordinate directions are provided. Please notice that + !! the cuboidal patch DOES NOT allow for the smearing of its + !! boundaries. + !! @param patch_id is the patch identifier + !! @param ib_markers_sf Array to track patch ids + subroutine s_ib_cuboid(patch_id, ib_markers_sf) + + integer, intent(in) :: patch_id + integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + + integer :: i, j, k !< Generic loop iterators + + ! Transferring the cuboid's centroid and length information + x_centroid = patch_ib(patch_id)%x_centroid + y_centroid = patch_ib(patch_id)%y_centroid + z_centroid = patch_ib(patch_id)%z_centroid + length_x = patch_ib(patch_id)%length_x + length_y = patch_ib(patch_id)%length_y + length_z = patch_ib(patch_id)%length_z + + ! Computing the beginning and the end x-, y- and z-coordinates of + ! the cuboid based on its centroid and lengths + x_boundary%beg = x_centroid - 0.5_wp*length_x + x_boundary%end = x_centroid + 0.5_wp*length_x + y_boundary%beg = y_centroid - 0.5_wp*length_y + y_boundary%end = y_centroid + 0.5_wp*length_y + z_boundary%beg = z_centroid - 0.5_wp*length_z + z_boundary%end = z_centroid + 0.5_wp*length_z + + ! Since the cuboidal patch does not allow for its boundaries to get + ! smoothed out, the pseudo volume fraction is set to 1 to make sure + ! that only the current patch contributes to the fluid state in the + ! cells that this patch covers. + eta = 1._wp + + ! Checking whether the cuboid covers a particular cell in the domain + ! and verifying whether the current patch has permission to write to + ! to that cell. If both queries check out, the primitive variables + ! of the current patch are assigned to this cell. + do k = 0, p + do j = 0, n + do i = 0, m + + if (grid_geometry == 3) then + call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) + else + cart_y = y_cc(j) + cart_z = z_cc(k) + end if + + if (x_boundary%beg <= x_cc(i) .and. & + x_boundary%end >= x_cc(i) .and. & + y_boundary%beg <= cart_y .and. & + y_boundary%end >= cart_y .and. & + z_boundary%beg <= cart_z .and. & + z_boundary%end >= cart_z) then + + ! Updating the patch identities bookkeeping variable + ib_markers_sf(i, j, k) = patch_id + end if + end do + end do + end do + + end subroutine s_ib_cuboid + + !> The cylindrical patch is a 3D geometry that may be used, + !! for example, in setting up a cylindrical solid boundary + !! confinement, like a blood vessel. The geometry of this + !! patch is well-defined when the centroid, the radius and + !! the length along the cylinder's axis, parallel to the x-, + !! y- or z-coordinate direction, are provided. Please note + !! that the cylindrical patch DOES allow for the smoothing + !! of its lateral boundary. + !! @param patch_id is the patch identifier + !! @param ib_markers_sf Array to track patch ids + !! @param ib True if this patch is an immersed boundary + subroutine s_ib_cylinder(patch_id, ib_markers_sf) + + integer, intent(in) :: patch_id + integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + + integer :: i, j, k !< Generic loop iterators + real(wp) :: radius + + ! Transferring the cylindrical patch's centroid, length, radius, + ! smoothing patch identity and smoothing coefficient information + + x_centroid = patch_ib(patch_id)%x_centroid + y_centroid = patch_ib(patch_id)%y_centroid + z_centroid = patch_ib(patch_id)%z_centroid + length_x = patch_ib(patch_id)%length_x + length_y = patch_ib(patch_id)%length_y + length_z = patch_ib(patch_id)%length_z + radius = patch_ib(patch_id)%radius + + ! Computing the beginning and the end x-, y- and z-coordinates of + ! the cylinder based on its centroid and lengths + x_boundary%beg = x_centroid - 0.5_wp*length_x + x_boundary%end = x_centroid + 0.5_wp*length_x + y_boundary%beg = y_centroid - 0.5_wp*length_y + y_boundary%end = y_centroid + 0.5_wp*length_y + z_boundary%beg = z_centroid - 0.5_wp*length_z + z_boundary%end = z_centroid + 0.5_wp*length_z + + ! Initializing the pseudo volume fraction value to 1. The value will + ! be modified as the patch is laid out on the grid, but only in the + ! case that smearing of the cylindrical patch's boundary is enabled. + eta = 1._wp + + ! Checking whether the cylinder covers a particular cell in the + ! domain and verifying whether the current patch has the permission + ! to write to that cell. If both queries check out, the primitive + ! variables of the current patch are assigned to this cell. + do k = 0, p + do j = 0, n + do i = 0, m + + if (grid_geometry == 3) then + call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) + else + cart_y = y_cc(j) + cart_z = z_cc(k) + end if + + if (((.not. f_is_default(length_x) .and. & + (cart_y - y_centroid)**2 & + + (cart_z - z_centroid)**2 <= radius**2 .and. & + x_boundary%beg <= x_cc(i) .and. & + x_boundary%end >= x_cc(i)) & + .or. & + (.not. f_is_default(length_y) .and. & + (x_cc(i) - x_centroid)**2 & + + (cart_z - z_centroid)**2 <= radius**2 .and. & + y_boundary%beg <= cart_y .and. & + y_boundary%end >= cart_y) & + .or. & + (.not. f_is_default(length_z) .and. & + (x_cc(i) - x_centroid)**2 & + + (cart_y - y_centroid)**2 <= radius**2 .and. & + z_boundary%beg <= cart_z .and. & + z_boundary%end >= cart_z))) then + + ! Updating the patch identities bookkeeping variable + ib_markers_sf(i, j, k) = patch_id + end if + end do + end do + end do + + end subroutine s_ib_cylinder + + !> The STL patch is a 2/3D geometry that is imported from an STL file. + !! @param patch_id is the patch identifier + !! @param ib_markers_sf Array to track patch ids + !! @param ib True if this patch is an immersed boundary + !! @param STL_levelset STL levelset + !! @param STL_levelset_norm STL levelset normals + subroutine s_ib_model(patch_id, ib_markers_sf, STL_levelset, STL_levelset_norm) + + integer, intent(in) :: patch_id + integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + + ! Variables for IBM+STL + type(levelset_field), optional, intent(inout) :: STL_levelset !< Levelset determined by models + type(levelset_norm_field), optional, intent(inout) :: STL_levelset_norm !< Levelset_norm determined by models + real(wp) :: normals(1:3) !< Boundary normal buffer + integer :: boundary_vertex_count, boundary_edge_count, total_vertices !< Boundary vertex + real(wp), allocatable, dimension(:, :, :) :: boundary_v !< Boundary vertex buffer + real(wp), allocatable, dimension(:, :) :: interpolated_boundary_v !< Interpolated vertex buffer + real(wp) :: distance !< Levelset distance buffer + logical :: interpolate !< Logical variable to determine whether or not the model should be interpolated + + integer :: i, j, k !< Generic loop iterators + + type(t_bbox) :: bbox, bbox_old + type(t_model) :: model + type(ic_model_parameters) :: params + + real(wp), dimension(1:3) :: point, model_center + + real(wp) :: grid_mm(1:3, 1:2) + + integer :: cell_num + integer :: ncells + + real(wp), dimension(1:4, 1:4) :: transform, transform_n + + print *, " * Reading model: "//trim(patch_ib(patch_id)%model_filepath) + + model = f_model_read(patch_ib(patch_id)%model_filepath) + params%scale(:) = patch_ib(patch_id)%model_scale(:) + params%translate(:) = patch_ib(patch_id)%model_translate(:) + params%rotate(:) = patch_ib(patch_id)%model_rotate(:) + params%spc = patch_ib(patch_id)%model_spc + params%threshold = patch_ib(patch_id)%model_threshold + + if (proc_rank == 0) then + print *, " * Transforming model." + end if + + ! Get the model center before transforming the model + bbox_old = f_create_bbox(model) + model_center(1:3) = (bbox_old%min(1:3) + bbox_old%max(1:3))/2._wp + + ! Compute the transform matrices for vertices and normals + transform = f_create_transform_matrix(params, model_center) + transform_n = f_create_transform_matrix(params) + + call s_transform_model(model, transform, transform_n) + + ! Recreate the bounding box after transformation + bbox = f_create_bbox(model) + + ! Show the number of vertices in the original STL model + if (proc_rank == 0) then + print *, ' * Number of input model vertices:', 3*model%ntrs + end if + + call f_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count) + + ! Check if the model needs interpolation + if (p > 0) then + call f_check_interpolation_3D(model, (/dx, dy, dz/), interpolate) + else + call f_check_interpolation_2D(boundary_v, boundary_edge_count, (/dx, dy, dz/), interpolate) + end if + + ! Show the number of edges and boundary edges in 2D STL models + if (proc_rank == 0 .and. p == 0) then + print *, ' * Number of 2D model boundary edges:', boundary_edge_count + end if + + ! Interpolate the STL model along the edges (2D) and on triangle facets (3D) + if (interpolate) then + if (proc_rank == 0) then + print *, ' * Interpolating STL vertices.' + end if + + if (p > 0) then + call f_interpolate_3D(model, (/dx, dy, dz/), interpolated_boundary_v, total_vertices) + else + call f_interpolate_2D(boundary_v, boundary_edge_count, (/dx, dy, dz/), interpolated_boundary_v, total_vertices) + end if + + if (proc_rank == 0) then + print *, ' * Total number of interpolated boundary vertices:', total_vertices + end if + end if + + if (proc_rank == 0) then + write (*, "(A, 3(2X, F20.10))") " > Model: Min:", bbox%min(1:3) + write (*, "(A, 3(2X, F20.10))") " > Cen:", (bbox%min(1:3) + bbox%max(1:3))/2._wp + write (*, "(A, 3(2X, F20.10))") " > Max:", bbox%max(1:3) + + !call s_model_write("__out__.stl", model) + !call s_model_write("__out__.obj", model) + + grid_mm(1, :) = (/minval(x_cc) - 0.e5_wp*dx, maxval(x_cc) + 0.e5_wp*dx/) + grid_mm(2, :) = (/minval(y_cc) - 0.e5_wp*dy, maxval(y_cc) + 0.e5_wp*dy/) + + if (p > 0) then + grid_mm(3, :) = (/minval(z_cc) - 0.e5_wp*dz, maxval(z_cc) + 0.e5_wp*dz/) + else + grid_mm(3, :) = (/0._wp, 0._wp/) + end if + + write (*, "(A, 3(2X, F20.10))") " > Domain: Min:", grid_mm(:, 1) + write (*, "(A, 3(2X, F20.10))") " > Cen:", (grid_mm(:, 1) + grid_mm(:, 2))/2._wp + write (*, "(A, 3(2X, F20.10))") " > Max:", grid_mm(:, 2) + end if + + ncells = (m + 1)*(n + 1)*(p + 1) + do i = 0, m; do j = 0, n; do k = 0, p + + cell_num = i*(n + 1)*(p + 1) + j*(p + 1) + (k + 1) + if (proc_rank == 0 .and. mod(cell_num, ncells/100) == 0) then + write (*, "(A, I3, A)", advance="no") & + char(13)//" * Generating grid: ", & + nint(100*real(cell_num)/ncells), "%" + end if + + point = (/x_cc(i), y_cc(j), 0._wp/) + if (p > 0) then + point(3) = z_cc(k) + end if + + if (grid_geometry == 3) then + point = f_convert_cyl_to_cart(point) + end if + + eta = f_model_is_inside(model, point, (/dx, dy, dz/), patch_ib(patch_id)%model_spc) + + ! Reading STL boundary vertices and compute the levelset and levelset_norm + if (eta > patch_ib(patch_id)%model_threshold) then + ib_markers_sf(i, j, k) = patch_id + end if + + ! 3D models + if (p > 0) then + + ! Get the boundary normals and shortest distance between the cell center and the model boundary + call f_distance_normals_3D(model, point, normals, distance) + + ! Get the shortest distance between the cell center and the interpolated model boundary + if (interpolate) then + STL_levelset%sf(i, j, k, patch_id) = f_interpolated_distance(interpolated_boundary_v, & + total_vertices, & + point) + else + STL_levelset%sf(i, j, k, patch_id) = distance + end if + + ! Correct the sign of the levelset + if (ib_markers_sf(i, j, k) > 0) then + STL_levelset%sf(i, j, k, patch_id) = -abs(STL_levelset%sf(i, j, k, patch_id)) + end if + + ! Correct the sign of the levelset_norm + if (ib_markers_sf(i, j, k) == 0) then + normals(1:3) = -normals(1:3) + end if + + ! Assign the levelset_norm + STL_levelset_norm%sf(i, j, k, patch_id, 1:3) = normals(1:3) + else + ! 2D models + if (interpolate) then + ! Get the shortest distance between the cell center and the model boundary + STL_levelset%sf(i, j, 0, patch_id) = f_interpolated_distance(interpolated_boundary_v, & + total_vertices, & + point) + else + ! Get the shortest distance between the cell center and the interpolated model boundary + STL_levelset%sf(i, j, 0, patch_id) = f_distance(boundary_v, & + boundary_edge_count, & + point) + end if + + ! Correct the sign of the levelset + if (ib_markers_sf(i, j, k) > 0) then + STL_levelset%sf(i, j, 0, patch_id) = -abs(STL_levelset%sf(i, j, 0, patch_id)) + end if + + ! Get the boundary normals + call f_normals(boundary_v, & + boundary_edge_count, & + point, & + normals) + + ! Correct the sign of the levelset_norm + if (ib_markers_sf(i, j, k) == 0) then + normals(1:3) = -normals(1:3) + end if + + ! Assign the levelset_norm + STL_levelset_norm%sf(i, j, k, patch_id, 1:3) = normals(1:3) + + end if + end do; end do; end do + + if (proc_rank == 0) then + print *, "" + print *, " * Cleaning up." + end if + + call s_model_free(model) + + end subroutine s_ib_model + + subroutine s_convert_cylindrical_to_cartesian_coord(cyl_y, cyl_z) + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: cyl_y, cyl_z + + cart_y = cyl_y*sin(cyl_z) + cart_z = cyl_y*cos(cyl_z) + + end subroutine s_convert_cylindrical_to_cartesian_coord + + pure function f_convert_cyl_to_cart(cyl) result(cart) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), dimension(1:3), intent(in) :: cyl + real(wp), dimension(1:3) :: cart + + cart = (/cyl(1), & + cyl(2)*sin(cyl(3)), & + cyl(2)*cos(cyl(3))/) + + end function f_convert_cyl_to_cart + + subroutine s_convert_cylindrical_to_spherical_coord(cyl_x, cyl_y) + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(IN) :: cyl_x, cyl_y + + sph_phi = atan(cyl_y/cyl_x) + + end subroutine s_convert_cylindrical_to_spherical_coord + + !> Archimedes spiral function + !! @param myth Angle + !! @param offset Thickness + !! @param a Starting position + pure elemental function f_r(myth, offset, a) + $:GPU_ROUTINE(parallelism='[seq]') + real(wp), intent(in) :: myth, offset, a + real(wp) :: b + real(wp) :: f_r + + !r(th) = a + b*th + + b = 2._wp*a/(2._wp*pi) + f_r = a + b*myth + offset + end function f_r + +end module m_ib_patches diff --git a/src/pre_process/m_model.fpp b/src/common/m_model.fpp similarity index 100% rename from src/pre_process/m_model.fpp rename to src/common/m_model.fpp diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 4332681f11..07463137f2 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -98,16 +98,6 @@ contains #ifdef MFC_MPI integer :: ierr !< Generic flag used to identify and report MPI errors -#endif - -#ifndef MFC_MPI - - ! Serial run only has 1 processor - num_procs = 1 - ! Local processor rank is 0 - proc_rank = 0 - -#else ! Initializing the MPI environment call MPI_INIT(ierr) @@ -123,7 +113,11 @@ contains ! Querying the rank of the local processor call MPI_COMM_RANK(MPI_COMM_WORLD, proc_rank, ierr) - +#else + ! Serial run only has 1 processor + num_procs = 1 + ! Local processor rank is 0 + proc_rank = 0 #endif end subroutine s_mpi_initialize @@ -168,27 +162,20 @@ contains end if !Additional variables pb and mv for non-polytropic qbmm -#ifdef MFC_PRE_PROCESS if (qbmm .and. .not. polytropic) then do i = 1, nb do j = 1, nnode +#ifdef MFC_PRE_PROCESS MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j)%sf => pb%sf(0:m, 0:n, 0:p, j, i) MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j + nb*nnode)%sf => mv%sf(0:m, 0:n, 0:p, j, i) - end do - end do - end if -#endif - -#ifdef MFC_SIMULATION - if (qbmm .and. .not. polytropic) then - do i = 1, nb - do j = 1, nnode +#elif defined (MFC_SIMULATION) MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j)%sf => pb_ts(1)%sf(0:m, 0:n, 0:p, j, i) MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j + nb*nnode)%sf => mv_ts(1)%sf(0:m, 0:n, 0:p, j, i) +#endif end do end do end if -#endif + ! Define global(g) and local(l) sizes for flow variables sizes_glb(1) = m_glb + 1; sizes_loc(1) = m + 1 if (n > 0) then @@ -485,6 +472,30 @@ contains end subroutine s_mpi_allreduce_sum + !> The following subroutine takes the input local variable + !! from all processors and reduces to the sum of all + !! values. The reduced variable is recorded back onto the + !! original local variable on each processor. + !! @param var_loc Some variable containing the local value which should be + !! reduced amongst all the processors in the communicator. + !! @param var_glb The globally reduced value + impure subroutine s_mpi_allreduce_integer_sum(var_loc, var_glb) + + integer, intent(in) :: var_loc + integer, intent(out) :: var_glb + +#ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors + + ! Performing the reduction procedure + call MPI_ALLREDUCE(var_loc, var_glb, 1, MPI_INTEGER, & + MPI_SUM, MPI_COMM_WORLD, ierr) +#else + var_glb = var_loc +#endif + + end subroutine s_mpi_allreduce_integer_sum + !> The following subroutine takes the input local variable !! from all processors and reduces to the minimum of all !! values. The reduced variable is recorded back onto the diff --git a/src/pre_process/m_checker.fpp b/src/pre_process/m_checker.fpp index 0001444584..f07e25dbfa 100644 --- a/src/pre_process/m_checker.fpp +++ b/src/pre_process/m_checker.fpp @@ -254,6 +254,10 @@ contains "Incompatible BC type for boundary condition patch "//trim(iStr)) end do - end subroutine + end subroutine s_check_bc + + impure subroutine s_check_moving_IBM + + end subroutine s_check_moving_IBM end module m_checker diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 0eb0cf6c41..ca379da0cc 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -538,6 +538,12 @@ contains patch_ib(i)%model_filepath(:) = dflt_char patch_ib(i)%model_spc = num_ray patch_ib(i)%model_threshold = ray_tracing_threshold + + ! Variables to handle moving imersed boundaries, defaulting to no movement + patch_ib(i)%moving_ibm = 0 + patch_ib(i)%vel(1) = 0._wp + patch_ib(i)%vel(2) = 0._wp + patch_ib(i)%vel(3) = 0._wp end do ! Fluids physical parameters diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_icpp_patches.fpp similarity index 64% rename from src/pre_process/m_patches.fpp rename to src/pre_process/m_icpp_patches.fpp index 071a34d996..0157be3767 100644 --- a/src/pre_process/m_patches.fpp +++ b/src/pre_process/m_icpp_patches.fpp @@ -9,7 +9,7 @@ #:include '3dHardcodedIC.fpp' #:include 'macros.fpp' -module m_patches +module m_icpp_patches use m_model ! Subroutine(s) related to STL files @@ -29,9 +29,11 @@ module m_patches use m_mpi_common + use m_ib_patches + implicit none - private; public :: s_apply_domain_patches + private; public :: s_apply_icpp_patches real(wp) :: x_centroid, y_centroid, z_centroid real(wp) :: length_x, length_y, length_z @@ -64,13 +66,10 @@ module m_patches contains - impure subroutine s_apply_domain_patches(patch_id_fp, q_prim_vf, ib_markers_sf, levelset, levelset_norm) + impure subroutine s_apply_icpp_patches(patch_id_fp, q_prim_vf) type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf integer, dimension(0:m, 0:m, 0:m), intent(inout) :: patch_id_fp - integer, dimension(:, :, :), intent(inout), optional :: ib_markers_sf - type(levelset_field), intent(inout), optional :: levelset !< Levelset determined by models - type(levelset_norm_field), intent(inout), optional :: levelset_norm !< Levelset_norm determined by models integer :: i @@ -87,55 +86,28 @@ contains !> @{ ! Spherical patch if (patch_icpp(i)%geometry == 8) then - call s_sphere(i, patch_id_fp, q_prim_vf) + call s_icpp_sphere(i, patch_id_fp, q_prim_vf) ! Cuboidal patch elseif (patch_icpp(i)%geometry == 9) then - call s_cuboid(i, patch_id_fp, q_prim_vf) + call s_icpp_cuboid(i, patch_id_fp, q_prim_vf) ! Cylindrical patch elseif (patch_icpp(i)%geometry == 10) then - call s_cylinder(i, patch_id_fp, q_prim_vf) + call s_icpp_cylinder(i, patch_id_fp, q_prim_vf) ! Swept plane patch elseif (patch_icpp(i)%geometry == 11) then - call s_sweep_plane(i, patch_id_fp, q_prim_vf) + call s_icpp_sweep_plane(i, patch_id_fp, q_prim_vf) ! Ellipsoidal patch elseif (patch_icpp(i)%geometry == 12) then - call s_ellipsoid(i, patch_id_fp, q_prim_vf) + call s_icpp_ellipsoid(i, patch_id_fp, q_prim_vf) ! Spherical harmonic patch elseif (patch_icpp(i)%geometry == 14) then - call s_spherical_harmonic(i, patch_id_fp, q_prim_vf) + call s_icpp_spherical_harmonic(i, patch_id_fp, q_prim_vf) ! 3D Modified circular patch elseif (patch_icpp(i)%geometry == 19) then - call s_3dvarcircle(i, patch_id_fp, q_prim_vf) + call s_icpp_3dvarcircle(i, patch_id_fp, q_prim_vf) ! 3D STL patch elseif (patch_icpp(i)%geometry == 21) then - call s_model(i, patch_id_fp, q_prim_vf) - end if - end do - !> @} - - !> IB Patches - !> @{ - ! Spherical patch - do i = 1, num_ibs - if (proc_rank == 0) then - print *, 'Processing 3D ib patch ', i - end if - - if (patch_ib(i)%geometry == 8) then - call s_sphere(i, ib_markers_sf, q_prim_vf, ib) - call s_sphere_levelset(i, levelset, levelset_norm) - elseif (patch_ib(i)%geometry == 9) then - call s_cuboid(i, ib_markers_sf, q_prim_vf, ib) - call s_cuboid_levelset(i, levelset, levelset_norm) - elseif (patch_ib(i)%geometry == 10) then - call s_cylinder(i, ib_markers_sf, q_prim_vf, ib) - call s_cylinder_levelset(i, levelset, levelset_norm) - elseif (patch_ib(i)%geometry == 11) then - call s_3D_airfoil(i, ib_markers_sf, q_prim_vf, ib) - call s_3D_airfoil_levelset(i, levelset, levelset_norm) - ! STL+IBM patch - elseif (patch_ib(i)%geometry == 12) then - call s_model(i, ib_markers_sf, q_prim_vf, ib, levelset, levelset_norm) + call s_icpp_model(i, patch_id_fp, q_prim_vf) end if end do !> @} @@ -153,61 +125,39 @@ contains !> @{ ! Circular patch if (patch_icpp(i)%geometry == 2) then - call s_circle(i, patch_id_fp, q_prim_vf) + call s_icpp_circle(i, patch_id_fp, q_prim_vf) ! Rectangular patch elseif (patch_icpp(i)%geometry == 3) then - call s_rectangle(i, patch_id_fp, q_prim_vf) + call s_icpp_rectangle(i, patch_id_fp, q_prim_vf) ! Swept line patch elseif (patch_icpp(i)%geometry == 4) then - call s_sweep_line(i, patch_id_fp, q_prim_vf) + call s_icpp_sweep_line(i, patch_id_fp, q_prim_vf) ! Elliptical patch elseif (patch_icpp(i)%geometry == 5) then - call s_ellipse(i, patch_id_fp, q_prim_vf) + call s_icpp_ellipse(i, patch_id_fp, q_prim_vf) ! Unimplemented patch (formerly isentropic vortex) elseif (patch_icpp(i)%geometry == 6) then call s_mpi_abort('This used to be the isentropic vortex patch, '// & 'which no longer exists. See Examples. Exiting.') ! Spherical Harmonic Patch elseif (patch_icpp(i)%geometry == 14) then - call s_spherical_harmonic(i, patch_id_fp, q_prim_vf) + call s_icpp_spherical_harmonic(i, patch_id_fp, q_prim_vf) ! Spiral patch elseif (patch_icpp(i)%geometry == 17) then - call s_spiral(i, patch_id_fp, q_prim_vf) + call s_icpp_spiral(i, patch_id_fp, q_prim_vf) ! Modified circular patch elseif (patch_icpp(i)%geometry == 18) then - call s_varcircle(i, patch_id_fp, q_prim_vf) + call s_icpp_varcircle(i, patch_id_fp, q_prim_vf) ! TaylorGreen vortex patch elseif (patch_icpp(i)%geometry == 20) then - call s_2D_TaylorGreen_vortex(i, patch_id_fp, q_prim_vf) + call s_icpp_2D_TaylorGreen_vortex(i, patch_id_fp, q_prim_vf) ! STL patch elseif (patch_icpp(i)%geometry == 21) then - call s_model(i, patch_id_fp, q_prim_vf) + call s_icpp_model(i, patch_id_fp, q_prim_vf) end if !> @} end do - !> IB Patches - !> @{ - do i = 1, num_ibs - if (proc_rank == 0) then - print *, 'Processing 2D ib patch ', i - end if - if (patch_ib(i)%geometry == 2) then - call s_circle(i, ib_markers_sf, q_prim_vf, ib) - call s_circle_levelset(i, levelset, levelset_norm) - elseif (patch_ib(i)%geometry == 3) then - call s_rectangle(i, ib_markers_sf, q_prim_vf, ib) - call s_rectangle_levelset(i, levelset, levelset_norm) - elseif (patch_ib(i)%geometry == 4) then - call s_airfoil(i, ib_markers_sf, q_prim_vf, ib) - call s_airfoil_levelset(i, levelset, levelset_norm) - ! STL+IBM patch - elseif (patch_ib(i)%geometry == 5) then - call s_model(i, ib_markers_sf, q_prim_vf, ib, levelset, levelset_norm) - end if - end do - !> @} - ! 1D Patch Geometries else @@ -219,16 +169,16 @@ contains ! Line segment patch if (patch_icpp(i)%geometry == 1) then - call s_line_segment(i, patch_id_fp, q_prim_vf) + call s_icpp_line_segment(i, patch_id_fp, q_prim_vf) ! 1d analytical elseif (patch_icpp(i)%geometry == 16) then - call s_1d_bubble_pulse(i, patch_id_fp, q_prim_vf) + call s_icpp_1d_bubble_pulse(i, patch_id_fp, q_prim_vf) end if end do end if - end subroutine s_apply_domain_patches + end subroutine s_apply_icpp_patches !> The line segment patch is a 1D geometry that may be used, !! for example, in creating a Riemann problem. The geometry @@ -239,7 +189,7 @@ contains !! @param patch_id patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_line_segment(patch_id, patch_id_fp, q_prim_vf) + subroutine s_icpp_line_segment(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp @@ -300,7 +250,7 @@ contains end do @:HardcodedDellacation() - end subroutine s_line_segment + end subroutine s_icpp_line_segment !> The spiral patch is a 2D geometry that may be used, The geometry !! of the patch is well-defined when its centroid and radius @@ -309,7 +259,7 @@ contains !! @param patch_id patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - impure subroutine s_spiral(patch_id, patch_id_fp, q_prim_vf) + impure subroutine s_icpp_spiral(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp @@ -370,7 +320,7 @@ contains end do @:HardcodedDellacation() - end subroutine s_spiral + end subroutine s_icpp_spiral !> The circular patch is a 2D geometry that may be used, for !! example, in creating a bubble or a droplet. The geometry @@ -380,13 +330,11 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - !! @param ib True if this patch is an immersed boundary - subroutine s_circle(patch_id, patch_id_fp, q_prim_vf, ib_flag) + subroutine s_icpp_circle(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib_flag real(wp) :: radius @@ -397,17 +345,11 @@ contains ! Transferring the circular patch's radius, centroid, smearing patch ! identity and smearing coefficient information - if (present(ib_flag)) then - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid - radius = patch_ib(patch_id)%radius - else - x_centroid = patch_icpp(patch_id)%x_centroid - y_centroid = patch_icpp(patch_id)%y_centroid - radius = patch_icpp(patch_id)%radius - smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id - smooth_coeff = patch_icpp(patch_id)%smooth_coeff - end if + x_centroid = patch_icpp(patch_id)%x_centroid + y_centroid = patch_icpp(patch_id)%y_centroid + radius = patch_icpp(patch_id)%radius + smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id + smooth_coeff = patch_icpp(patch_id)%smooth_coeff ! Initializing the pseudo volume fraction value to 1. The value will ! be modified as the patch is laid out on the grid, but only in the @@ -422,7 +364,7 @@ contains do j = 0, n do i = 0, m - if (.not. present(ib_flag) .and. patch_icpp(patch_id)%smoothen) then + if (patch_icpp(patch_id)%smoothen) then eta = tanh(smooth_coeff/min(dx, dy)* & (sqrt((x_cc(i) - x_centroid)**2 & @@ -431,376 +373,35 @@ contains end if - if (present(ib_flag) .and. ((x_cc(i) - x_centroid)**2 & - + (y_cc(j) - y_centroid)**2 <= radius**2)) & + if (((x_cc(i) - x_centroid)**2 & + + (y_cc(j) - y_centroid)**2 <= radius**2 & + .and. & + patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) & + .or. & + patch_id_fp(i, j, 0) == smooth_patch_id) & then - patch_id_fp(i, j, 0) = patch_id - else - if (((x_cc(i) - x_centroid)**2 & - + (y_cc(j) - y_centroid)**2 <= radius**2 & - .and. & - patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) & - .or. & - (.not. present(ib_flag) .and. patch_id_fp(i, j, 0) == smooth_patch_id)) & - then - - call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - eta, q_prim_vf, patch_id_fp) - - @:analytical() - if (patch_icpp(patch_id)%hcid /= dflt_int) then - @:Hardcoded2D() - end if - - end if - end if - end do - end do - @:HardcodedDellacation() - - end subroutine s_circle - - !! @param patch_id is the patch identifier - !! @param patch_id_fp Array to track patch ids - !! @param q_prim_vf Array of primitive variables - !! @param ib True if this patch is an immersed boundary - subroutine s_airfoil(patch_id, patch_id_fp, q_prim_vf, ib_flag) - - integer, intent(in) :: patch_id - integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp - type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib_flag - - real(wp) :: x0, y0, f, x_act, y_act, ca_in, pa, ma, ta, theta - real(wp) :: xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c - integer :: i, j, k - integer :: Np1, Np2 - - if (.not. present(ib_flag)) return - x0 = patch_ib(patch_id)%x_centroid - y0 = patch_ib(patch_id)%y_centroid - ca_in = patch_ib(patch_id)%c - pa = patch_ib(patch_id)%p - ma = patch_ib(patch_id)%m - ta = patch_ib(patch_id)%t - theta = pi*patch_ib(patch_id)%theta/180._wp - - Np1 = int((pa*ca_in/dx)*20) - Np2 = int(((ca_in - pa*ca_in)/dx)*20) - Np = Np1 + Np2 + 1 - - allocate (airfoil_grid_u(1:Np)) - allocate (airfoil_grid_l(1:Np)) - - airfoil_grid_u(1)%x = x0 - airfoil_grid_u(1)%y = y0 - - airfoil_grid_l(1)%x = x0 - airfoil_grid_l(1)%y = y0 - - eta = 1._wp - - do i = 1, Np1 + Np2 - 1 - if (i <= Np1) then - xc = x0 + i*(pa*ca_in/Np1) - xa = (xc - x0)/ca_in - yc = (ma/pa**2)*(2*pa*xa - xa**2) - dycdxc = (2*ma/pa**2)*(pa - xa) - else - xc = x0 + pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2) - xa = (xc - x0)/ca_in - yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) - dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) - end if - - yt = (5._wp*ta)*(0.2969_wp*xa**0.5_wp - 0.126_wp*xa - 0.3516_wp*xa**2._wp + 0.2843_wp*xa**3 - 0.1015_wp*xa**4) - sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp - cos_c = 1/(1 + dycdxc**2)**0.5_wp - - xu = xa - yt*sin_c - yu = yc + yt*cos_c - - xl = xa + yt*sin_c - yl = yc - yt*cos_c - - xu = xu*ca_in + x0 - yu = yu*ca_in + y0 - - xl = xl*ca_in + x0 - yl = yl*ca_in + y0 - - airfoil_grid_u(i + 1)%x = xu - airfoil_grid_u(i + 1)%y = yu - - airfoil_grid_l(i + 1)%x = xl - airfoil_grid_l(i + 1)%y = yl - - end do - - airfoil_grid_u(Np)%x = x0 + ca_in - airfoil_grid_u(Np)%y = y0 - - airfoil_grid_l(Np)%x = x0 + ca_in - airfoil_grid_l(Np)%y = y0 - - do j = 0, n - do i = 0, m - - if (.not. f_is_default(patch_ib(patch_id)%theta)) then - x_act = (x_cc(i) - x0)*cos(theta) - (y_cc(j) - y0)*sin(theta) + x0 - y_act = (x_cc(i) - x0)*sin(theta) + (y_cc(j) - y0)*cos(theta) + y0 - else - x_act = x_cc(i) - y_act = y_cc(j) - end if + call s_assign_patch_primitive_variables(patch_id, i, j, 0, & + eta, q_prim_vf, patch_id_fp) - if (x_act >= x0 .and. x_act <= x0 + ca_in) then - xa = (x_act - x0)/ca_in - if (xa <= pa) then - yc = (ma/pa**2)*(2*pa*xa - xa**2) - dycdxc = (2*ma/pa**2)*(pa - xa) - else - yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) - dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) + @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded2D() end if - if (y_act >= y0) then - k = 1 - do while (airfoil_grid_u(k)%x < x_act) - k = k + 1 - end do - if (f_approx_equal(airfoil_grid_u(k)%x, x_act)) then - if (y_act <= airfoil_grid_u(k)%y) then - !!IB - !call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - !eta, q_prim_vf, patch_id_fp) - patch_id_fp(i, j, 0) = patch_id - end if - else - f = (airfoil_grid_u(k)%x - x_act)/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x) - if (y_act <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then - !!IB - !call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - !eta, q_prim_vf, patch_id_fp) - patch_id_fp(i, j, 0) = patch_id - end if - end if - else - k = 1 - do while (airfoil_grid_l(k)%x < x_act) - k = k + 1 - end do - if (f_approx_equal(airfoil_grid_l(k)%x, x_act)) then - if (y_act >= airfoil_grid_l(k)%y) then - !!IB - !call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - !eta, q_prim_vf, patch_id_fp) - patch_id_fp(i, j, 0) = patch_id - end if - else - f = (airfoil_grid_l(k)%x - x_act)/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x) - if (y_act >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then - !!IB - !call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - !eta, q_prim_vf, patch_id_fp) - patch_id_fp(i, j, 0) = patch_id - end if - end if - end if end if end do end do + @:HardcodedDellacation() - if (.not. f_is_default(patch_ib(patch_id)%theta)) then - do i = 1, Np - airfoil_grid_l(i)%x = (airfoil_grid_l(i)%x - x0)*cos(theta) + (airfoil_grid_l(i)%y - y0)*sin(theta) + x0 - airfoil_grid_l(i)%y = -1._wp*(airfoil_grid_l(i)%x - x0)*sin(theta) + (airfoil_grid_l(i)%y - y0)*cos(theta) + y0 - - airfoil_grid_u(i)%x = (airfoil_grid_u(i)%x - x0)*cos(theta) + (airfoil_grid_u(i)%y - y0)*sin(theta) + x0 - airfoil_grid_u(i)%y = -1._wp*(airfoil_grid_u(i)%x - x0)*sin(theta) + (airfoil_grid_u(i)%y - y0)*cos(theta) + y0 - end do - end if - - end subroutine s_airfoil - - !! @param patch_id is the patch identifier - !! @param patch_id_fp Array to track patch ids - !! @param q_prim_vf Array of primitive variables - !! @param ib True if this patch is an immersed boundary - subroutine s_3D_airfoil(patch_id, patch_id_fp, q_prim_vf, ib_flag) - - integer, intent(in) :: patch_id - integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp - type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib_flag - - real(wp) :: x0, y0, z0, lz, z_max, z_min, f, x_act, y_act, ca_in, pa, ma, ta, theta, xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c - integer :: i, j, k, l - integer :: Np1, Np2 - - if (.not. present(ib_flag)) return - x0 = patch_ib(patch_id)%x_centroid - y0 = patch_ib(patch_id)%y_centroid - z0 = patch_ib(patch_id)%z_centroid - lz = patch_ib(patch_id)%length_z - ca_in = patch_ib(patch_id)%c - pa = patch_ib(patch_id)%p - ma = patch_ib(patch_id)%m - ta = patch_ib(patch_id)%t - theta = pi*patch_ib(patch_id)%theta/180._wp - - Np1 = int((pa*ca_in/dx)*20) - Np2 = int(((ca_in - pa*ca_in)/dx)*20) - Np = Np1 + Np2 + 1 - - allocate (airfoil_grid_u(1:Np)) - allocate (airfoil_grid_l(1:Np)) - - airfoil_grid_u(1)%x = x0 - airfoil_grid_u(1)%y = y0 - - airfoil_grid_l(1)%x = x0 - airfoil_grid_l(1)%y = y0 - - z_max = z0 + lz/2 - z_min = z0 - lz/2 - - eta = 1._wp - - do i = 1, Np1 + Np2 - 1 - if (i <= Np1) then - xc = x0 + i*(pa*ca_in/Np1) - xa = (xc - x0)/ca_in - yc = (ma/pa**2)*(2*pa*xa - xa**2) - dycdxc = (2*ma/pa**2)*(pa - xa) - else - xc = x0 + pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2) - xa = (xc - x0)/ca_in - yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) - dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) - end if - - yt = (5._wp*ta)*(0.2969_wp*xa**0.5_wp - 0.126_wp*xa - 0.3516_wp*xa**2._wp + 0.2843_wp*xa**3 - 0.1015_wp*xa**4) - sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp - cos_c = 1/(1 + dycdxc**2)**0.5_wp - - xu = xa - yt*sin_c - yu = yc + yt*cos_c - - xl = xa + yt*sin_c - yl = yc - yt*cos_c - - xu = xu*ca_in + x0 - yu = yu*ca_in + y0 - - xl = xl*ca_in + x0 - yl = yl*ca_in + y0 - - airfoil_grid_u(i + 1)%x = xu - airfoil_grid_u(i + 1)%y = yu - - airfoil_grid_l(i + 1)%x = xl - airfoil_grid_l(i + 1)%y = yl - - end do - - airfoil_grid_u(Np)%x = x0 + ca_in - airfoil_grid_u(Np)%y = y0 - - airfoil_grid_l(Np)%x = x0 + ca_in - airfoil_grid_l(Np)%y = y0 - - do l = 0, p - if (z_cc(l) >= z_min .and. z_cc(l) <= z_max) then - do j = 0, n - do i = 0, m - - if (.not. f_is_default(patch_ib(patch_id)%theta)) then - x_act = (x_cc(i) - x0)*cos(theta) - (y_cc(j) - y0)*sin(theta) + x0 - y_act = (x_cc(i) - x0)*sin(theta) + (y_cc(j) - y0)*cos(theta) + y0 - else - x_act = x_cc(i) - y_act = y_cc(j) - end if - - if (x_act >= x0 .and. x_act <= x0 + ca_in) then - xa = (x_act - x0)/ca_in - if (xa <= pa) then - yc = (ma/pa**2)*(2*pa*xa - xa**2) - dycdxc = (2*ma/pa**2)*(pa - xa) - else - yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) - dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) - end if - if (y_act >= y0) then - k = 1 - do while (airfoil_grid_u(k)%x < x_act) - k = k + 1 - end do - if (f_approx_equal(airfoil_grid_u(k)%x, x_act)) then - if (y_act <= airfoil_grid_u(k)%y) then - !!IB - !call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - !eta, q_prim_vf, patch_id_fp) - patch_id_fp(i, j, l) = patch_id - end if - else - f = (airfoil_grid_u(k)%x - x_act)/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x) - if (y_act <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then - !!IB - !call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - !eta, q_prim_vf, patch_id_fp) - patch_id_fp(i, j, l) = patch_id - end if - end if - else - k = 1 - do while (airfoil_grid_l(k)%x < x_act) - k = k + 1 - end do - if (f_approx_equal(airfoil_grid_l(k)%x, x_act)) then - if (y_act >= airfoil_grid_l(k)%y) then - !!IB - !call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - !eta, q_prim_vf, patch_id_fp) - patch_id_fp(i, j, l) = patch_id - end if - else - f = (airfoil_grid_l(k)%x - x_act)/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x) - - if (y_act >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then - !!IB - !call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - !eta, q_prim_vf, patch_id_fp) - patch_id_fp(i, j, l) = patch_id - end if - end if - end if - end if - end do - end do - end if - end do - - if (.not. f_is_default(patch_ib(patch_id)%theta)) then - do i = 1, Np - airfoil_grid_l(i)%x = (airfoil_grid_l(i)%x - x0)*cos(theta) + (airfoil_grid_l(i)%y - y0)*sin(theta) + x0 - airfoil_grid_l(i)%y = -1._wp*(airfoil_grid_l(i)%x - x0)*sin(theta) + (airfoil_grid_l(i)%y - y0)*cos(theta) + y0 - - airfoil_grid_u(i)%x = (airfoil_grid_u(i)%x - x0)*cos(theta) + (airfoil_grid_u(i)%y - y0)*sin(theta) + x0 - airfoil_grid_u(i)%y = -1._wp*(airfoil_grid_u(i)%x - x0)*sin(theta) + (airfoil_grid_u(i)%y - y0)*cos(theta) + y0 - end do - end if - - end subroutine s_3D_airfoil + end subroutine s_icpp_circle !> The varcircle patch is a 2D geometry that may be used !! . It generatres an annulus !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_varcircle(patch_id, patch_id_fp, q_prim_vf) + subroutine s_icpp_varcircle(patch_id, patch_id_fp, q_prim_vf) ! Patch identifier integer, intent(in) :: patch_id @@ -859,12 +460,12 @@ contains end do @:HardcodedDellacation() - end subroutine s_varcircle + end subroutine s_icpp_varcircle !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_3dvarcircle(patch_id, patch_id_fp, q_prim_vf) + subroutine s_icpp_3dvarcircle(patch_id, patch_id_fp, q_prim_vf) ! Patch identifier integer, intent(in) :: patch_id @@ -929,7 +530,7 @@ contains end do @:HardcodedDellacation() - end subroutine s_3dvarcircle + end subroutine s_icpp_3dvarcircle !> The elliptical patch is a 2D geometry. The geometry of !! the patch is well-defined when its centroid and radii @@ -938,7 +539,7 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_ellipse(patch_id, patch_id_fp, q_prim_vf) + subroutine s_icpp_ellipse(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp @@ -1001,7 +602,7 @@ contains end do @:HardcodedDellacation() - end subroutine s_ellipse + end subroutine s_icpp_ellipse !> The ellipsoidal patch is a 3D geometry. The geometry of !! the patch is well-defined when its centroid and radii @@ -1010,7 +611,7 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_ellipsoid(patch_id, patch_id_fp, q_prim_vf) + subroutine s_icpp_ellipsoid(patch_id, patch_id_fp, q_prim_vf) ! Patch identifier integer, intent(in) :: patch_id @@ -1088,7 +689,7 @@ contains end do @:HardcodedDellacation() - end subroutine s_ellipsoid + end subroutine s_icpp_ellipsoid !> The rectangular patch is a 2D geometry that may be used, !! for example, in creating a solid boundary, or pre-/post- @@ -1101,13 +702,11 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - !! @param ib True if this patch is an immersed boundary - subroutine s_rectangle(patch_id, patch_id_fp, q_prim_vf, ib_flag) + subroutine s_icpp_rectangle(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib_flag !< True if this patch is an immersed boundary integer :: i, j, k !< generic loop iterators real(wp) :: pi_inf, gamma, lit_gamma !< Equation of state parameters @@ -1119,17 +718,10 @@ contains lit_gamma = (1._wp + gamma)/gamma ! Transferring the rectangle's centroid and length information - if (present(ib_flag)) then - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid - length_x = patch_ib(patch_id)%length_x - length_y = patch_ib(patch_id)%length_y - else - x_centroid = patch_icpp(patch_id)%x_centroid - y_centroid = patch_icpp(patch_id)%y_centroid - length_x = patch_icpp(patch_id)%length_x - length_y = patch_icpp(patch_id)%length_y - end if + x_centroid = patch_icpp(patch_id)%x_centroid + y_centroid = patch_icpp(patch_id)%y_centroid + length_x = patch_icpp(patch_id)%length_x + length_y = patch_icpp(patch_id)%length_y ! Computing the beginning and the end x- and y-coordinates of the ! rectangle based on its centroid and lengths @@ -1154,40 +746,34 @@ contains x_boundary%end >= x_cc(i) .and. & y_boundary%beg <= y_cc(j) .and. & y_boundary%end >= y_cc(j)) then - if (present(ib_flag)) then - ! Updating the patch identities bookkeeping variable - patch_id_fp(i, j, 0) = patch_id - else - if (patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) & - then - - call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - eta, q_prim_vf, patch_id_fp) - - @:analytical() + if (patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) & + then - if (patch_icpp(patch_id)%hcid /= dflt_int) then - @:Hardcoded2D() - end if + call s_assign_patch_primitive_variables(patch_id, i, j, 0, & + eta, q_prim_vf, patch_id_fp) - if ((q_prim_vf(1)%sf(i, j, 0) < 1.e-10) .and. (model_eqns == 4)) then - !zero density, reassign according to Tait EOS - q_prim_vf(1)%sf(i, j, 0) = & - (((q_prim_vf(E_idx)%sf(i, j, 0) + pi_inf)/(pref + pi_inf))**(1._wp/lit_gamma))* & - rhoref*(1._wp - q_prim_vf(alf_idx)%sf(i, j, 0)) - end if + @:analytical() - ! Updating the patch identities bookkeeping variable - if (1._wp - eta < sgm_eps) patch_id_fp(i, j, 0) = patch_id + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded2D() + end if + if ((q_prim_vf(1)%sf(i, j, 0) < 1.e-10) .and. (model_eqns == 4)) then + !zero density, reassign according to Tait EOS + q_prim_vf(1)%sf(i, j, 0) = & + (((q_prim_vf(E_idx)%sf(i, j, 0) + pi_inf)/(pref + pi_inf))**(1._wp/lit_gamma))* & + rhoref*(1._wp - q_prim_vf(alf_idx)%sf(i, j, 0)) end if + + ! Updating the patch identities bookkeeping variable + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, 0) = patch_id end if end if end do end do @:HardcodedDellacation() - end subroutine s_rectangle + end subroutine s_icpp_rectangle !> The swept line patch is a 2D geometry that may be used, !! for example, in creating a solid boundary, or pre-/post- @@ -1199,7 +785,7 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_sweep_line(patch_id, patch_id_fp, q_prim_vf) + subroutine s_icpp_sweep_line(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp @@ -1261,7 +847,7 @@ contains end do @:HardcodedDellacation() - end subroutine s_sweep_line + end subroutine s_icpp_sweep_line !> The Taylor Green vortex is 2D decaying vortex that may be used, !! for example, to verify the effects of viscous attenuation. @@ -1270,7 +856,7 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_2D_TaylorGreen_Vortex(patch_id, patch_id_fp, q_prim_vf) + subroutine s_icpp_2D_TaylorGreen_Vortex(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp @@ -1343,12 +929,12 @@ contains end do @:HardcodedDellacation() - end subroutine s_2D_TaylorGreen_Vortex + end subroutine s_icpp_2D_TaylorGreen_Vortex !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_1d_bubble_pulse(patch_id, patch_id_fp, q_prim_vf) + subroutine s_icpp_1d_bubble_pulse(patch_id, patch_id_fp, q_prim_vf) ! Description: This patch assigns the primitive variables as analytical ! functions such that the code can be verified. @@ -1404,14 +990,14 @@ contains end do @:HardcodedDellacation() - end subroutine s_1D_bubble_pulse + end subroutine s_icpp_1D_bubble_pulse !> This patch generates the shape of the spherical harmonics !! as a perturbation to a perfect sphere !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_spherical_harmonic(patch_id, patch_id_fp, q_prim_vf) + subroutine s_icpp_spherical_harmonic(patch_id, patch_id_fp, q_prim_vf) integer, intent(IN) :: patch_id integer, intent(INOUT), dimension(0:m, 0:n, 0:p) :: patch_id_fp @@ -1548,7 +1134,7 @@ contains end do end if - end subroutine s_spherical_harmonic + end subroutine s_icpp_spherical_harmonic !> The spherical patch is a 3D geometry that may be used, !! for example, in creating a bubble or a droplet. The patch @@ -1558,13 +1144,11 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - !! @param ib True if this patch is an immersed boundary - subroutine s_sphere(patch_id, patch_id_fp, q_prim_vf, ib_flag) + subroutine s_icpp_sphere(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib_flag !< True if this patch is an immersed boundary ! Generic loop iterators integer :: i, j, k @@ -1577,19 +1161,12 @@ contains ! Transferring spherical patch's radius, centroid, smoothing patch ! identity and smoothing coefficient information - if (present(ib_flag)) then - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid - z_centroid = patch_ib(patch_id)%z_centroid - radius = patch_ib(patch_id)%radius - else - x_centroid = patch_icpp(patch_id)%x_centroid - y_centroid = patch_icpp(patch_id)%y_centroid - z_centroid = patch_icpp(patch_id)%z_centroid - radius = patch_icpp(patch_id)%radius - smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id - smooth_coeff = patch_icpp(patch_id)%smooth_coeff - end if + x_centroid = patch_icpp(patch_id)%x_centroid + y_centroid = patch_icpp(patch_id)%y_centroid + z_centroid = patch_icpp(patch_id)%z_centroid + radius = patch_icpp(patch_id)%radius + smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id + smooth_coeff = patch_icpp(patch_id)%smooth_coeff ! Initializing the pseudo volume fraction value to 1. The value will ! be modified as the patch is laid out on the grid, but only in the @@ -1611,7 +1188,7 @@ contains cart_z = z_cc(k) end if - if (.not. present(ib_flag) .and. patch_icpp(patch_id)%smoothen) then + if (patch_icpp(patch_id)%smoothen) then eta = tanh(smooth_coeff/min(dx, dy, dz)* & (sqrt((x_cc(i) - x_centroid)**2 & + (cart_y - y_centroid)**2 & @@ -1619,36 +1196,27 @@ contains - radius))*(-0.5_wp) + 0.5_wp end if - if (present(ib_flag)) then - ! Updating the patch identities bookkeeping variable - if (((x_cc(i) - x_centroid)**2 & - + (cart_y - y_centroid)**2 & - + (cart_z - z_centroid)**2 <= radius**2)) then - patch_id_fp(i, j, k) = patch_id - end if - else - if ((((x_cc(i) - x_centroid)**2 & - + (cart_y - y_centroid)**2 & - + (cart_z - z_centroid)**2 <= radius**2) .and. & - patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. & - patch_id_fp(i, j, k) == smooth_patch_id) then - - call s_assign_patch_primitive_variables(patch_id, i, j, k, & - eta, q_prim_vf, patch_id_fp) + if ((((x_cc(i) - x_centroid)**2 & + + (cart_y - y_centroid)**2 & + + (cart_z - z_centroid)**2 <= radius**2) .and. & + patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. & + patch_id_fp(i, j, k) == smooth_patch_id) then - @:analytical() - if (patch_icpp(patch_id)%hcid /= dflt_int) then - @:Hardcoded3D() - end if + call s_assign_patch_primitive_variables(patch_id, i, j, k, & + eta, q_prim_vf, patch_id_fp) + @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded3D() end if + end if end do end do end do @:HardcodedDellacation() - end subroutine s_sphere + end subroutine s_icpp_sphere !> The cuboidal patch is a 3D geometry that may be used, for !! example, in creating a solid boundary, or pre-/post-shock @@ -1661,10 +1229,9 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_cuboid(patch_id, patch_id_fp, q_prim_vf, ib_flag) + subroutine s_icpp_cuboid(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id - logical, optional, intent(in) :: ib_flag integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf @@ -1673,21 +1240,12 @@ contains @:Hardcoded3DVariables() ! Transferring the cuboid's centroid and length information - if (present(ib_flag)) then - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid - z_centroid = patch_ib(patch_id)%z_centroid - length_x = patch_ib(patch_id)%length_x - length_y = patch_ib(patch_id)%length_y - length_z = patch_ib(patch_id)%length_z - else - x_centroid = patch_icpp(patch_id)%x_centroid - y_centroid = patch_icpp(patch_id)%y_centroid - z_centroid = patch_icpp(patch_id)%z_centroid - length_x = patch_icpp(patch_id)%length_x - length_y = patch_icpp(patch_id)%length_y - length_z = patch_icpp(patch_id)%length_z - end if + x_centroid = patch_icpp(patch_id)%x_centroid + y_centroid = patch_icpp(patch_id)%y_centroid + z_centroid = patch_icpp(patch_id)%z_centroid + length_x = patch_icpp(patch_id)%length_x + length_y = patch_icpp(patch_id)%length_y + length_z = patch_icpp(patch_id)%length_z ! Computing the beginning and the end x-, y- and z-coordinates of ! the cuboid based on its centroid and lengths @@ -1726,24 +1284,19 @@ contains z_boundary%beg <= cart_z .and. & z_boundary%end >= cart_z) then - if (present(ib_flag)) then - ! Updating the patch identities bookkeeping variable - patch_id_fp(i, j, k) = patch_id - else - if (patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) then + if (patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) then - call s_assign_patch_primitive_variables(patch_id, i, j, k, & - eta, q_prim_vf, patch_id_fp) + call s_assign_patch_primitive_variables(patch_id, i, j, k, & + eta, q_prim_vf, patch_id_fp) - @:analytical() - if (patch_icpp(patch_id)%hcid /= dflt_int) then - @:Hardcoded3D() - end if + @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded3D() + end if - ! Updating the patch identities bookkeeping variable - if (1._wp - eta < sgm_eps) patch_id_fp(i, j, k) = patch_id + ! Updating the patch identities bookkeeping variable + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, k) = patch_id - end if end if end if end do @@ -1751,7 +1304,7 @@ contains end do @:HardcodedDellacation() - end subroutine s_cuboid + end subroutine s_icpp_cuboid !> The cylindrical patch is a 3D geometry that may be used, !! for example, in setting up a cylindrical solid boundary @@ -1764,13 +1317,11 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - !! @param ib True if this patch is an immersed boundary - subroutine s_cylinder(patch_id, patch_id_fp, q_prim_vf, ib_flag) + subroutine s_icpp_cylinder(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib_flag !< True if this patch is an immersed boundary integer :: i, j, k !< Generic loop iterators real(wp) :: radius @@ -1779,26 +1330,15 @@ contains ! Transferring the cylindrical patch's centroid, length, radius, ! smoothing patch identity and smoothing coefficient information - - if (present(ib_flag)) then - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid - z_centroid = patch_ib(patch_id)%z_centroid - length_x = patch_ib(patch_id)%length_x - length_y = patch_ib(patch_id)%length_y - length_z = patch_ib(patch_id)%length_z - radius = patch_ib(patch_id)%radius - else - x_centroid = patch_icpp(patch_id)%x_centroid - y_centroid = patch_icpp(patch_id)%y_centroid - z_centroid = patch_icpp(patch_id)%z_centroid - length_x = patch_icpp(patch_id)%length_x - length_y = patch_icpp(patch_id)%length_y - length_z = patch_icpp(patch_id)%length_z - radius = patch_icpp(patch_id)%radius - smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id - smooth_coeff = patch_icpp(patch_id)%smooth_coeff - end if + x_centroid = patch_icpp(patch_id)%x_centroid + y_centroid = patch_icpp(patch_id)%y_centroid + z_centroid = patch_icpp(patch_id)%z_centroid + length_x = patch_icpp(patch_id)%length_x + length_y = patch_icpp(patch_id)%length_y + length_z = patch_icpp(patch_id)%length_z + radius = patch_icpp(patch_id)%radius + smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id + smooth_coeff = patch_icpp(patch_id)%smooth_coeff ! Computing the beginning and the end x-, y- and z-coordinates of ! the cylinder based on its centroid and lengths @@ -1829,7 +1369,7 @@ contains cart_z = z_cc(k) end if - if (.not. present(ib_flag) .and. patch_icpp(patch_id)%smoothen) then + if (patch_icpp(patch_id)%smoothen) then if (.not. f_is_default(length_x)) then eta = tanh(smooth_coeff/min(dy, dz)* & (sqrt((cart_y - y_centroid)**2 & @@ -1848,68 +1388,43 @@ contains end if end if - if (present(ib_flag)) then - if (((.not. f_is_default(length_x) .and. & - (cart_y - y_centroid)**2 & - + (cart_z - z_centroid)**2 <= radius**2 .and. & - x_boundary%beg <= x_cc(i) .and. & - x_boundary%end >= x_cc(i)) & - .or. & - (.not. f_is_default(length_y) .and. & - (x_cc(i) - x_centroid)**2 & - + (cart_z - z_centroid)**2 <= radius**2 .and. & - y_boundary%beg <= cart_y .and. & - y_boundary%end >= cart_y) & - .or. & - (.not. f_is_default(length_z) .and. & - (x_cc(i) - x_centroid)**2 & - + (cart_y - y_centroid)**2 <= radius**2 .and. & - z_boundary%beg <= cart_z .and. & - z_boundary%end >= cart_z))) then - - ! Updating the patch identities bookkeeping variable - patch_id_fp(i, j, k) = patch_id - end if + if (((.not. f_is_default(length_x) .and. & + (cart_y - y_centroid)**2 & + + (cart_z - z_centroid)**2 <= radius**2 .and. & + x_boundary%beg <= x_cc(i) .and. & + x_boundary%end >= x_cc(i)) & + .or. & + (.not. f_is_default(length_y) .and. & + (x_cc(i) - x_centroid)**2 & + + (cart_z - z_centroid)**2 <= radius**2 .and. & + y_boundary%beg <= cart_y .and. & + y_boundary%end >= cart_y) & + .or. & + (.not. f_is_default(length_z) .and. & + (x_cc(i) - x_centroid)**2 & + + (cart_y - y_centroid)**2 <= radius**2 .and. & + z_boundary%beg <= cart_z .and. & + z_boundary%end >= cart_z) .and. & + patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. & + patch_id_fp(i, j, k) == smooth_patch_id) then - else - if (((.not. f_is_default(length_x) .and. & - (cart_y - y_centroid)**2 & - + (cart_z - z_centroid)**2 <= radius**2 .and. & - x_boundary%beg <= x_cc(i) .and. & - x_boundary%end >= x_cc(i)) & - .or. & - (.not. f_is_default(length_y) .and. & - (x_cc(i) - x_centroid)**2 & - + (cart_z - z_centroid)**2 <= radius**2 .and. & - y_boundary%beg <= cart_y .and. & - y_boundary%end >= cart_y) & - .or. & - (.not. f_is_default(length_z) .and. & - (x_cc(i) - x_centroid)**2 & - + (cart_y - y_centroid)**2 <= radius**2 .and. & - z_boundary%beg <= cart_z .and. & - z_boundary%end >= cart_z) .and. & - patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. & - patch_id_fp(i, j, k) == smooth_patch_id) then - - call s_assign_patch_primitive_variables(patch_id, i, j, k, & - eta, q_prim_vf, patch_id_fp) - - @:analytical() - if (patch_icpp(patch_id)%hcid /= dflt_int) then - @:Hardcoded3D() - end if + call s_assign_patch_primitive_variables(patch_id, i, j, k, & + eta, q_prim_vf, patch_id_fp) - ! Updating the patch identities bookkeeping variable - if (1._wp - eta < sgm_eps) patch_id_fp(i, j, k) = patch_id + @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded3D() end if + + ! Updating the patch identities bookkeeping variable + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, k) = patch_id end if end do end do end do @:HardcodedDellacation() - end subroutine s_cylinder + end subroutine s_icpp_cylinder !> The swept plane patch is a 3D geometry that may be used, !! for example, in creating a solid boundary, or pre-/post- @@ -1921,7 +1436,7 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Primitive variables - subroutine s_sweep_plane(patch_id, patch_id_fp, q_prim_vf) + subroutine s_icpp_sweep_plane(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp @@ -1997,25 +1512,21 @@ contains end do @:HardcodedDellacation() - end subroutine s_sweep_plane + end subroutine s_icpp_sweep_plane !> The STL patch is a 2/3D geometry that is imported from an STL file. !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Primitive variables - !! @param ib True if this patch is an immersed boundary !! @param STL_levelset STL levelset !! @param STL_levelset_norm STL levelset normals - subroutine s_model(patch_id, patch_id_fp, q_prim_vf, ib_flag, STL_levelset, STL_levelset_norm) + subroutine s_icpp_model(patch_id, patch_id_fp, q_prim_vf) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf ! Variables for IBM+STL - type(levelset_field), optional, intent(inout) :: STL_levelset !< Levelset determined by models - type(levelset_norm_field), optional, intent(inout) :: STL_levelset_norm !< Levelset_norm determined by models - logical, optional, intent(in) :: ib_flag !< True if this patch is an immersed boundary real(wp) :: normals(1:3) !< Boundary normal buffer integer :: boundary_vertex_count, boundary_edge_count, total_vertices !< Boundary vertex real(wp), allocatable, dimension(:, :, :) :: boundary_v !< Boundary vertex buffer @@ -2038,27 +1549,16 @@ contains real(wp), dimension(1:4, 1:4) :: transform, transform_n - if (present(ib_flag) .and. proc_rank == 0) then - print *, " * Reading model: "//trim(patch_ib(patch_id)%model_filepath) - else if (proc_rank == 0) then + if (proc_rank == 0) then print *, " * Reading model: "//trim(patch_icpp(patch_id)%model_filepath) end if - if (present(ib_flag)) then - model = f_model_read(patch_ib(patch_id)%model_filepath) - params%scale(:) = patch_ib(patch_id)%model_scale(:) - params%translate(:) = patch_ib(patch_id)%model_translate(:) - params%rotate(:) = patch_ib(patch_id)%model_rotate(:) - params%spc = patch_ib(patch_id)%model_spc - params%threshold = patch_ib(patch_id)%model_threshold - else - model = f_model_read(patch_icpp(patch_id)%model_filepath) - params%scale(:) = patch_icpp(patch_id)%model_scale(:) - params%translate(:) = patch_icpp(patch_id)%model_translate(:) - params%rotate(:) = patch_icpp(patch_id)%model_rotate(:) - params%spc = patch_icpp(patch_id)%model_spc - params%threshold = patch_icpp(patch_id)%model_threshold - end if + model = f_model_read(patch_icpp(patch_id)%model_filepath) + params%scale(:) = patch_icpp(patch_id)%model_scale(:) + params%translate(:) = patch_icpp(patch_id)%model_translate(:) + params%rotate(:) = patch_icpp(patch_id)%model_rotate(:) + params%spc = patch_icpp(patch_id)%model_spc + params%threshold = patch_icpp(patch_id)%model_threshold if (proc_rank == 0) then print *, " * Transforming model." @@ -2154,98 +1654,25 @@ contains point = f_convert_cyl_to_cart(point) end if - if (present(ib_flag)) then - eta = f_model_is_inside(model, point, (/dx, dy, dz/), patch_ib(patch_id)%model_spc) - else - eta = f_model_is_inside(model, point, (/dx, dy, dz/), patch_icpp(patch_id)%model_spc) - end if - - if (present(ib_flag)) then - ! Reading STL boundary vertices and compute the levelset and levelset_norm - if (eta > patch_ib(patch_id)%model_threshold) then - patch_id_fp(i, j, k) = patch_id - end if - - ! 3D models - if (p > 0) then - - ! Get the boundary normals and shortest distance between the cell center and the model boundary - call f_distance_normals_3D(model, point, normals, distance) - - ! Get the shortest distance between the cell center and the interpolated model boundary - if (interpolate) then - STL_levelset%sf(i, j, k, patch_id) = f_interpolated_distance(interpolated_boundary_v, & - total_vertices, & - point) - else - STL_levelset%sf(i, j, k, patch_id) = distance - end if - - ! Correct the sign of the levelset - if (patch_id_fp(i, j, k) > 0) then - STL_levelset%sf(i, j, k, patch_id) = -abs(STL_levelset%sf(i, j, k, patch_id)) - end if - - ! Correct the sign of the levelset_norm - if (patch_id_fp(i, j, k) == 0) then - normals(1:3) = -normals(1:3) - end if - - ! Assign the levelset_norm - STL_levelset_norm%sf(i, j, k, patch_id, 1:3) = normals(1:3) - else - ! 2D models - if (interpolate) then - ! Get the shortest distance between the cell center and the model boundary - STL_levelset%sf(i, j, 0, patch_id) = f_interpolated_distance(interpolated_boundary_v, & - total_vertices, & - point) - else - ! Get the shortest distance between the cell center and the interpolated model boundary - STL_levelset%sf(i, j, 0, patch_id) = f_distance(boundary_v, & - boundary_edge_count, & - point) - end if - - ! Correct the sign of the levelset - if (patch_id_fp(i, j, k) > 0) then - STL_levelset%sf(i, j, 0, patch_id) = -abs(STL_levelset%sf(i, j, 0, patch_id)) - end if - - ! Get the boundary normals - call f_normals(boundary_v, & - boundary_edge_count, & - point, & - normals) - - ! Correct the sign of the levelset_norm - if (patch_id_fp(i, j, k) == 0) then - normals(1:3) = -normals(1:3) - end if - - ! Assign the levelset_norm - STL_levelset_norm%sf(i, j, k, patch_id, 1:3) = normals(1:3) + eta = f_model_is_inside(model, point, (/dx, dy, dz/), patch_icpp(patch_id)%model_spc) + if (patch_icpp(patch_id)%smoothen) then + if (eta > patch_icpp(patch_id)%model_threshold) then + eta = 1._wp end if else - if (patch_icpp(patch_id)%smoothen) then - if (eta > patch_icpp(patch_id)%model_threshold) then - eta = 1._wp - end if + if (eta > patch_icpp(patch_id)%model_threshold) then + eta = 1._wp else - if (eta > patch_icpp(patch_id)%model_threshold) then - eta = 1._wp - else - eta = 0._wp - end if + eta = 0._wp end if - call s_assign_patch_primitive_variables(patch_id, i, j, k, & - eta, q_prim_vf, patch_id_fp) - - ! Note: Should probably use *eta* to compute primitive variables - ! if defining them analytically. - @:analytical() end if + call s_assign_patch_primitive_variables(patch_id, i, j, k, & + eta, q_prim_vf, patch_id_fp) + + ! Note: Should probably use *eta* to compute primitive variables + ! if defining them analytically. + @:analytical() end do; end do; end do if (proc_rank == 0) then @@ -2255,7 +1682,7 @@ contains call s_model_free(model) - end subroutine s_model + end subroutine s_icpp_model subroutine s_convert_cylindrical_to_cartesian_coord(cyl_y, cyl_z) $:GPU_ROUTINE(parallelism='[seq]') @@ -2305,4 +1732,4 @@ contains f_r = a + b*myth + offset end function f_r -end module m_patches +end module m_icpp_patches diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp index 53fc9811a2..bf3963ee74 100644 --- a/src/pre_process/m_initial_condition.fpp +++ b/src/pre_process/m_initial_condition.fpp @@ -27,7 +27,9 @@ module m_initial_condition use m_variables_conversion ! Subroutines to change the state variables from ! one form to another - use m_patches + use m_ib_patches + + use m_icpp_patches use m_assign_variables @@ -188,10 +190,9 @@ contains end if if (ib) then - call s_apply_domain_patches(patch_id_fp, q_prim_vf, ib_markers%sf, levelset, levelset_norm) - else - call s_apply_domain_patches(patch_id_fp, q_prim_vf) + call s_apply_ib_patches(ib_markers%sf, levelset, levelset_norm) end if + call s_apply_icpp_patches(patch_id_fp, q_prim_vf) if (num_bc_patches > 0) call s_apply_boundary_patches(q_prim_vf, bc_type) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 9e8b9f70fe..85b3d18bd5 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -28,7 +28,9 @@ module m_start_up use m_compile_specific !< Compile-specific procedures - use m_patches + use m_ib_patches + + use m_icpp_patches use m_assign_variables diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 4d01608120..6beb433f8e 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -22,6 +22,10 @@ module m_ibm use m_constants + use m_compute_levelset + + use m_ib_patches + implicit none private :: s_compute_image_points, & @@ -34,6 +38,7 @@ module m_ibm s_ibm_correct_state, & s_finalize_ibm_module + integer, allocatable, dimension(:, :, :) :: patch_id_fp type(integer_field), public :: ib_markers type(levelset_field), public :: levelset type(levelset_norm_field), public :: levelset_norm @@ -47,6 +52,8 @@ module m_ibm integer :: num_inner_gps !< Number of ghost points $:GPU_DECLARE(create='[gp_layers,num_gps,num_inner_gps]') + logical :: moving_immersed_boundary_flag + contains !> Allocates memory for the variables in the IBM module @@ -81,6 +88,24 @@ contains impure subroutine s_ibm_setup() integer :: i, j, k + integer :: max_num_gps, max_num_inner_gps + + moving_immersed_boundary_flag = .false. + do i = 1, num_ibs + if (patch_ib(i)%moving_ibm .ne. 0) then + moving_immersed_boundary_flag = .true. + exit + end if + end do + + ! $:GPU_UPDATE(device='[patch_ib]') + ! $:GPU_UPDATE(host='[patch_ib]') + ! print *, "Moving IBM ", patch_ib(1)%moving_ibm + ! print *, "X Velocity ", patch_ib(1)%vel(1) + ! print *, "Y Velocity ", patch_ib(1)%vel(2) + + ! Allocating the patch identities bookkeeping variable + allocate (patch_id_fp(0:m, 0:n, 0:p)) $:GPU_UPDATE(device='[ib_markers%sf]') $:GPU_UPDATE(device='[levelset%sf]') @@ -91,11 +116,14 @@ contains $:GPU_UPDATE(host='[ib_markers%sf]') + ! find the number of ghost points and set them to be the maximum total across ranks call s_find_num_ghost_points(num_gps, num_inner_gps) + call s_mpi_allreduce_integer_sum(num_gps, max_num_gps) + call s_mpi_allreduce_integer_sum(num_inner_gps, max_num_inner_gps) $:GPU_UPDATE(device='[num_gps, num_inner_gps]') - @:ALLOCATE(ghost_points(1:num_gps)) - @:ALLOCATE(inner_points(1:num_inner_gps)) + @:ALLOCATE(ghost_points(1:int(max_num_gps * 1.2))) + @:ALLOCATE(inner_points(1:int(max_num_inner_gps * 1.2))) $:GPU_ENTER_DATA(copyin='[ghost_points,inner_points]') @@ -239,7 +267,15 @@ contains vel_norm_IP = sum(vel_IP*norm)*norm vel_g = vel_IP - vel_norm_IP else - vel_g = 0._wp + if (patch_ib(patch_id)%moving_ibm == 0) then + ! we know the object is not moving if moving_ibm is 0 (false) + vel_g = 0._wp + else + do q = 1, 3 + ! if mibm is 1 or 2, then the boundary may be moving + vel_g(q) = patch_ib(patch_id)%vel(q) + end do + end if end if ! Set momentum @@ -391,6 +427,7 @@ contains end if if (f_approx_equal(norm(dim), 0._wp)) then + ! if the ghost point is almost equal to a cell location, we set it equal and continue ghost_points_in(q)%ip_grid(dim) = ghost_points_in(q)%loc(dim) else if (norm(dim) > 0) then @@ -405,6 +442,8 @@ contains .or. temp_loc > s_cc(index + 1))) index = index + dir if (index < -buff_size .or. index > bound) then + print *, q, index, bound, buff_size + print *, "temp_loc=", temp_loc, " s_cc(index)=", s_cc(index), " s_cc(index+1)=", s_cc(index + 1) print *, "Increase buff_size further in m_helper_basic (currently set to a minimum of 10)" error stop "Increase buff_size" end if @@ -471,7 +510,7 @@ contains end subroutine s_find_num_ghost_points !> Function that finds the ghost points - pure subroutine s_find_ghost_points(ghost_points_in, inner_points_in) + subroutine s_find_ghost_points(ghost_points_in, inner_points_in) type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points_in type(ghost_point), dimension(num_inner_gps), intent(INOUT) :: inner_points_in @@ -489,6 +528,7 @@ contains do i = 0, m do j = 0, n if (p == 0) then + ! 2D if (ib_markers%sf(i, j, 0) /= 0) then subsection_2D = ib_markers%sf( & i - gp_layers:i + gp_layers, & @@ -498,6 +538,7 @@ contains patch_id = ib_markers%sf(i, j, 0) ghost_points_in(count)%ib_patch_id = & patch_id + ghost_points_in(count)%slip = patch_ib(patch_id)%slip ! ghost_points(count)%rank = proc_rank @@ -530,6 +571,7 @@ contains end if end if else + ! 3D do k = 0, p if (ib_markers%sf(i, j, k) /= 0) then subsection_3D = ib_markers%sf( & @@ -860,6 +902,66 @@ contains end subroutine s_interpolate_image_point + !> Subroutine the updates the moving imersed boundary positions via Euler's method + impure subroutine s_propagate_mib(patch_id) + + integer, intent(in) :: patch_id + integer :: i + + ! start by using euler's method naiively, but eventually incorporate more sophistocation + if (patch_ib(patch_id)%moving_ibm .eq. 1) then + ! this continues with euler's method, which is obviously not that great and we need to add acceleration + do i = 1, 3 + patch_ib(patch_id)%vel(i) = patch_ib(patch_id)%vel(i) + 0.0*dt ! TODO :: ADD EXTERNAL FORCES HERE + end do + + patch_ib(patch_id)%x_centroid = patch_ib(patch_id)%x_centroid + patch_ib(patch_id)%vel(1)*dt + patch_ib(patch_id)%y_centroid = patch_ib(patch_id)%y_centroid + patch_ib(patch_id)%vel(2)*dt + patch_ib(patch_id)%z_centroid = patch_ib(patch_id)%z_centroid + patch_ib(patch_id)%vel(3)*dt + end if + + end subroutine s_propagate_mib + + !> Resets the current indexes of immersed boundaries and replaces them after updating + !> the position of each moving immersed boundary + impure subroutine s_update_mib(num_ibs, levelset, levelset_norm) + + integer, intent(in) :: num_ibs + type(levelset_field), intent(inout) :: levelset + type(levelset_norm_field), intent(inout) :: levelset_norm + + integer :: i + + ! Clears the existing immersed boundary indices + ib_markers%sf = 0 + + do i = 1, num_ibs + if (patch_ib(i)%moving_ibm .ne. 0) then + call s_propagate_mib(i) ! TODO :: THIS IS DONE TERRIBLY WITH EULER METHOD + end if + end do + $:GPU_UPDATE(device='[patch_ib]') + + ! recompute the new ib_patch locations and broadcast them. + call s_apply_ib_patches(ib_markers%sf(0:m, 0:n, 0:p), levelset, levelset_norm) + call s_populate_ib_buffers() ! transmitts the new IB markers via MPI + $:GPU_UPDATE(device='[ib_markers%sf]') + + ! recalculate the ghost point locations and coefficients + call s_find_num_ghost_points(num_gps, num_inner_gps) + $:GPU_UPDATE(device='[num_gps, num_inner_gps]') + + call s_find_ghost_points(ghost_points, inner_points) + ! $:GPU_UPDATE(device='[ghost_points, inner_points]') + + call s_compute_image_points(ghost_points, levelset, levelset_norm) + ! $:GPU_UPDATE(device='[ghost_points]') + + call s_compute_interpolation_coeffs(ghost_points) + $:GPU_UPDATE(device='[ghost_points]') + + end subroutine s_update_mib + !> Subroutine to deallocate memory reserved for the IBM module impure subroutine s_finalize_ibm_module() diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 93e0126ff8..617afb7c98 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -203,7 +203,8 @@ contains do i = 1, num_ibs #:for VAR in [ 'radius', 'length_x', 'length_y', & - & 'x_centroid', 'y_centroid', 'c', 'm', 'p', 't', 'theta', 'slip' ] + & 'x_centroid', 'y_centroid', 'c', 'm', 'p', 't', 'theta', 'slip', & + 'moving_ibm', 'vel',] call MPI_BCAST(patch_ib(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) #:endfor call MPI_BCAST(patch_ib(i)%geometry, 2, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) diff --git a/src/simulation/p_main.fpp b/src/simulation/p_main.fpp index 29cb3b8281..92f92f357b 100644 --- a/src/simulation/p_main.fpp +++ b/src/simulation/p_main.fpp @@ -37,6 +37,7 @@ program p_main call nvtxStartRange("INIT") + !Initialize MPI call nvtxStartRange("INIT-MPI") call s_initialize_mpi_domain() @@ -74,6 +75,10 @@ program p_main ! Time-stepping Loop do + if (moving_immersed_boundary_flag) then + call s_update_mib(num_ibs, levelset, levelset_norm) + end if + if (cfl_dt) then if (mytime >= t_stop) then call s_save_performance_metrics(time_avg, time_final, io_time_avg, & diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index 836811eefa..fb30a21f61 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -7,6 +7,7 @@ from .state import ARG from .run import case_dicts + QPVF_IDX_VARS = { 'alpha_rho': 'contxb', 'vel' : 'momxb', 'pres': 'E_idx', 'alpha': 'advxb', 'tau_e': 'stress_idx%beg', 'Y': 'chemxb', @@ -234,7 +235,7 @@ def __get_sim_fpp(self, print: bool) -> str: igr_pres_lim = 1 if self.params.get("igr_pres_lim", 'F') == 'T' else 0 # Throw error if wenoz_q is required but not set - return f"""\ + out = f"""\ #:set MFC_CASE_OPTIMIZATION = {ARG("case_optimization")} #:set recon_type = {recon_type} #:set weno_order = {weno_order} @@ -262,10 +263,11 @@ def __get_sim_fpp(self, print: bool) -> str: #:set viscous = {viscous} """ - return """\ -! This file is purposefully empty. It is only important for builds that make use -! of --case-optimization. -""" + else: + out = "" + + # We need to also include the pre_processing includes so that common subroutines have access to the @:analytical function + return out + f"\n{self.__get_pre_fpp(print)}" def get_fpp(self, target, print = True) -> str: def _prepend() -> str: @@ -277,7 +279,7 @@ def _default(_) -> str: return "! This file is purposefully empty." result = { - "pre_process" : self.__get_pre_fpp, + "pre_process" : self.__get_pre_fpp, "simulation" : self.__get_sim_fpp, }.get(build.get_target(target).name, _default)(print) diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 8378d3044d..c2d5efbeff 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -112,9 +112,13 @@ def analytic(self): for real_attr, ty in [("geometry", ParamType.INT), ("radius", ParamType.REAL), ("theta", ParamType.REAL), ("slip", ParamType.LOG), ("c", ParamType.REAL), ("p", ParamType.REAL), - ("t", ParamType.REAL), ("m", ParamType.REAL)]: + ("t", ParamType.REAL), ("m", ParamType.REAL), + ("moving_ibm", ParamType.INT)]: PRE_PROCESS[f"patch_ib({ib_id})%{real_attr}"] = ty + for vel_id in range(1, 4): + PRE_PROCESS[f"patch_ib({ib_id})%vel({vel_id})"] = ParamType.REAL + for cmp_id, cmp in enumerate(["x", "y", "z"]): cmp_id += 1 PRE_PROCESS[f'patch_ib({ib_id})%{cmp}_centroid'] = ParamType.REAL @@ -339,9 +343,13 @@ def analytic(self): for real_attr, ty in [("geometry", ParamType.INT), ("radius", ParamType.REAL), ("theta", ParamType.REAL), ("slip", ParamType.LOG), ("c", ParamType.REAL), ("p", ParamType.REAL), - ("t", ParamType.REAL), ("m", ParamType.REAL)]: + ("t", ParamType.REAL), ("m", ParamType.REAL), + ("moving_ibm", ParamType.INT)]: SIMULATION[f"patch_ib({ib_id})%{real_attr}"] = ty + for vel_id in range(1, 4): + SIMULATION[f"patch_ib({ib_id})%vel({vel_id})"] = ParamType.REAL + for cmp_id, cmp in enumerate(["x", "y", "z"]): cmp_id += 1 SIMULATION[f'patch_ib({ib_id})%{cmp}_centroid'] = ParamType.REAL diff --git a/toolchain/mfc/run/run.py b/toolchain/mfc/run/run.py index fb7d528ff9..687e8034ea 100644 --- a/toolchain/mfc/run/run.py +++ b/toolchain/mfc/run/run.py @@ -136,7 +136,7 @@ def __execute_job_script(qsystem: queues.QueueSystem): def run(targets = None, case = None): targets = get_targets(list(REQUIRED_TARGETS) + (targets or ARG("targets"))) case = case or input.load(ARG("input"), ARG("--")) - + build(targets) cons.print("[bold]Run[/bold]")