Resolve mibm error temp #1004
Closed
Resolve mibm error temp #1004
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
…ndaries with moving imersed boundary variables
…r what is applied each loop. Wrote a seaparate function and working on passing everything in
…ib patch definitions. This is a precommit before I start separating the two files
…ry patches from the m_ib_patches.fpp file
…ve compiling issues
…s crashing. Intermittent commit
…FDEF. There is likely a more-elegant way to do this in the future that involves unifying the defintion of dx in pre_process and simulation
Contributor
PR Reviewer Guide 🔍Here are some key observations to aid the review process:
|
![qodo-merge-pro[bot]](https://avatars.githubusercontent.com/in/484649?s=60&v=4){kind=small}
Contributor
There was a problem hiding this comment.
High-level Suggestion
The moving immersed boundary feature is inefficient because it recomputes all boundary data structures from scratch at every time step. This should be optimized to only update data for cells affected by the boundary's movement. [High-level, importance: 9]
Solution Walkthrough:
Before:
subroutine s_update_mib(...)
// Clear all existing IB markers on the entire grid
ib_markers%sf = 0
// Update position of all moving boundaries
for i = 1 to num_ibs:
call s_propagate_mib(i)
// Re-apply all IB patches to the entire grid from scratch
call s_apply_ib_patches(ib_markers%sf, ...)
// Recalculate all ghost points and their properties from scratch
call s_find_num_ghost_points(...)
call s_find_ghost_points(...)
call s_compute_image_points(...)
call s_compute_interpolation_coeffs(...)
end subroutine
After:
subroutine s_update_mib(...)
// For each moving boundary, identify the old and new regions
old_region = get_region(boundary_at_t)
new_region = get_region(boundary_at_t_plus_dt)
affected_region = union(old_region, new_region)
// Update IB markers only in the affected region
update_ib_markers_locally(affected_region)
// Identify ghost points that are added, removed, or need updates
// in the vicinity of the affected region.
affected_ghost_points = find_affected_ghost_points(affected_region)
// Recompute properties only for the affected ghost points
for each gp in affected_ghost_points:
recompute_image_point(gp)
recompute_interpolation_coeffs(gp)
end subroutine
Comment on lines
+376
to
392
| 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 |
Contributor
There was a problem hiding this comment.
Suggestion: Restore the missing logic to update the patch_id_fp array within the s_icpp_circle subroutine to ensure correct patch identification. [possible issue, importance: 8]
Suggested change
| 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 | |
| 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 | |
| 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 | |
| ! Updating the patch identities bookkeeping variable | |
| if (1._wp - eta < sgm_eps) patch_id_fp(i, j, 0) = patch_id | |
| end if |
Comment on lines
+1659
to
+1671
| 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) |
Contributor
There was a problem hiding this comment.
Suggestion: Correct the logic for handling smoothed and non-smoothed STL patches to ensure the eta value is calculated correctly based on the smoothen flag. [possible issue, importance: 8]
Suggested change
| 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) | |
| if (.not. patch_icpp(patch_id)%smoothen) then | |
| if (eta > patch_icpp(patch_id)%model_threshold) then | |
| eta = 1._wp | |
| else | |
| eta = 0._wp | |
| end if | |
| end if | |
| call s_assign_patch_primitive_variables(patch_id, i, j, k, & | |
| eta, q_prim_vf, patch_id_fp) |
Comment on lines
+911
to
+921
| ! 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 |
Contributor
There was a problem hiding this comment.
Suggestion: The velocity update logic is currently a no-op as it adds zero. Replace the placeholder 0.0 with the correct acceleration term or remove the line to reflect that velocity is constant until acceleration is implemented. [possible issue, importance: 5]
Suggested change
| ! 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 | |
| ! 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 | |
| ! TODO :: ADD EXTERNAL FORCES HERE to update velocity, e.g.: | |
| ! do i = 1, 3 | |
| ! patch_ib(patch_id)%vel(i) = patch_ib(patch_id)%vel(i) + acceleration(i)*dt | |
| ! 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 |
| bound = p + buff_size - 1 | ||
| end if | ||
|
|
||
| if (f_approx_equal(norm(dim), 0._wp)) then |
Contributor
There was a problem hiding this comment.
Suggestion: Prevent a potential out-of-bounds array access by adding an index < bound check to the do while loop condition, ensuring s_cc(index + 1) is accessed safely. [possible issue, importance: 8]
Member
|
haven't looked at all details yet but sniff test seems like this is quite good. i have a separate question i'll ask on slack |
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1004 +/- ##
==========================================
- Coverage 41.91% 41.74% -0.18%
==========================================
Files 69 70 +1
Lines 19783 20002 +219
Branches 2473 2483 +10
==========================================
+ Hits 8293 8349 +56
- Misses 9952 10121 +169
+ Partials 1538 1532 -6 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
Member
Author
|
This was accidentally submitted as a debug branch. A new PR has been created that has all of these changes. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
User description
Description
This PR adds a 1st-order moving immersed boundary method, which will need to be iterated upon. This implimented will require some additional work via adding rotations, higher-order time stepping, and performance optimizations. But it puts in the framework for moving immersed boundaries that can be built upon. I would like to move forward with merging, and then going back to optimize the code, since the current functionality completely captures the scope of a new feature. The code still performs optimally when the new feature is disabled.
This PR also refactors the patch frame work to separate the notions of ICPP and IB patches, which allows IB patches to be moved into a common directory.
Type of change
Please delete options that are not relevant.
Scope
If you cannot check the above box, please split your PR into multiple PRs that each have a common goal.
How Has This Been Tested?
Please describe the tests that you ran to verify your changes.
Provide instructions so we can reproduce.
Please also list any relevant details for your test configuration
Test Configuration:
I tested this locally on GNU compilers with and without MPI enabled. I also tested on Wingtip with NVHPC compilers with and without MPI and on GPUs.
Checklist
docs/)examples/that demonstrate my new feature performing as expected.They run to completion and demonstrate "interesting physics"
./mfc.sh formatbefore committing my codeIf your code changes any code source files (anything in
src/simulation)To make sure the code is performing as expected on GPU devices, I have:
nvtxranges so that they can be identified in profiles./mfc.sh run XXXX --gpu -t simulation --nsys, and have attached the output file (.nsys-rep) and plain text results to this PR./mfc.sh run XXXX --gpu -t simulation --rsys --hip-trace, and have attached the output file and plain text results to this PR.PR Type
Enhancement
Description
• Implements first-order moving immersed boundary method (MIBM) with velocity-based boundary conditions
• Refactors patch framework to separate ICPP and IB patches into dedicated modules (
m_icpp_patchesandm_ib_patches)• Adds moving boundary support with
moving_ibmflag and velocity parameters in patch data structures• Integrates moving boundary updates in main simulation time-stepping loop via
s_update_mibfunction• Enhances MPI support for broadcasting moving boundary parameters across processes
• Updates build configuration to exclude IB-specific modules from post-process builds
• Adds analytical function support in case generation toolchain
Diagram Walkthrough
File Walkthrough
13 files
m_icpp_patches.fpp
Refactor patches module to separate ICPP and IB functionalitysrc/pre_process/m_icpp_patches.fpp
• Renamed module from
m_patchestom_icpp_patchesand refactored toseparate ICPP and IB patch handling
• Removed IB patch processing
logic and parameters from all geometry subroutines
• Added
s_icpp_prefix to all geometry subroutines (sphere, circle, rectangle, etc.)
•
Simplified subroutine signatures by removing optional IB-related
parameters
m_ibm.fpp
Implement moving immersed boundary method functionalitysrc/simulation/m_ibm.fpp
• Added moving immersed boundary (MIBM) support with velocity-based
boundary conditions
• Implemented
s_propagate_mibands_update_mibfunctions for moving boundary updates
• Added
moving_immersed_boundary_flagto track if any boundaries are moving•
Enhanced ghost point allocation with buffer scaling for dynamic
boundaries
m_mpi_common.fpp
Add MPI integer reduction and fix compilation structuresrc/common/m_mpi_common.fpp
• Added
s_mpi_allreduce_integer_sumfunction for integer reductionoperations
• Fixed conditional compilation structure for MPI
initialization
• Consolidated QBMM variable handling across
pre-process and simulation modules
m_initial_condition.fpp
Update initial condition to use separated patch modulessrc/pre_process/m_initial_condition.fpp
• Updated to use separate
m_ib_patchesandm_icpp_patchesmodules•
Modified patch application to call both IB and ICPP patch functions
separately
m_mpi_proxy.fpp
Add MPI support for moving boundary parameterssrc/simulation/m_mpi_proxy.fpp
• Added MPI broadcasting for new IB patch parameters:
moving_ibmandvelp_main.fpp
Integrate moving boundary updates in main simulation loopsrc/simulation/p_main.fpp
• Added moving immersed boundary update call in main time-stepping
loop
• Integrated
s_update_mibfunction whenmoving_immersed_boundary_flagis activem_global_parameters.fpp
Initialize moving boundary parameters in global setupsrc/pre_process/m_global_parameters.fpp
• Added initialization of moving immersed boundary parameters
(
moving_ibm,vel) with default valuesm_derived_types.fpp
Add moving boundary data structures to IB patch typesrc/common/m_derived_types.fpp
• Added
moving_ibminteger flag andvelvelocity array toib_patch_parameterstypem_checker.fpp
Add placeholder for moving IBM validationsrc/pre_process/m_checker.fpp
• Added placeholder
s_check_moving_IBMsubroutine for futurevalidation
m_start_up.fpp
Update startup module imports for separated patchessrc/pre_process/m_start_up.fpp
• Updated module imports to use separated
m_ib_patchesandm_icpp_patchescase.py
Enhance case generation to support analytical functionstoolchain/mfc/case.py
• Modified FPP generation to include pre-processing includes in
simulation builds
• Added support for common subroutines to access
@:analyticalfunctioncase_dicts.py
Add moving boundary parameters to case dictionariestoolchain/mfc/run/case_dicts.py
• Added
moving_ibmparameter andvelvelocity components to IB patchparameter dictionaries
• Extended parameter support for both
pre-process and simulation modules
m_ib_patches.fpp
New immersed boundary patches module implementationsrc/common/m_ib_patches.fpp
• Creates a new module
m_ib_patchesdedicated to immersed boundary(IB) patch handling
• Implements geometry-specific IB patch
subroutines for circles, rectangles, spheres, cuboids, cylinders,
airfoils, and STL models
• Adds coordinate conversion utilities for
cylindrical to cartesian and spherical coordinates
• Includes levelset
computation functionality for various geometric shapes
2 files
case.fpp
Minor formatting change in case include filesrc/common/include/case.fpp
• Added empty line in analytical macro definition
run.py
Minor formatting change in run moduletoolchain/mfc/run/run.py
• Added trailing whitespace in run function
2 files
CMakeLists.txt
Update build configuration for separated modulesCMakeLists.txt
• Added exclusion of
m_compute_levelset.fppandm_ib_patches.fppfrompost_process builds
• Added
stdc++library linking for SILO buildssettings.json
Disable Fortran language server in VS Code settings.vscode/settings.json
• Disabled fortls language server by setting
fortran.fortls.disabledto true
6 files