Invicid two-way fluid-structure interaction #1075
There are no files selected for viewing
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
| Original file line number | Diff line number | Diff line change | ||||
|---|---|---|---|---|---|---|
| @@ -0,0 +1,105 @@ | ||||||
| import json | ||||||
| import math | ||||||
|
|
||||||
| Mu = 1.84e-05 | ||||||
|
There was a problem hiding this comment. Suggestion: Unused variable: Severity Level: Critical 🚨
Suggested change
Why it matters? ⭐The variable Mu is defined but not referenced anywhere in the file after the PR. It's dead code and should be removed to avoid confusion; this is a harmless cleanup that doesn't change behavior. Prompt for AI Agent 🤖This is a comment left during a code review.
**Path:** examples/2D_mibm_cylinder_in_cross_flow/case.py
**Line:** 4:4
**Comment:**
*Possible Bug: Unused variable: `Mu` is assigned but never referenced; this is dead code and can be removed to avoid misleading readers about intended viscous computations.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise. |
||||||
| gam_a = 1.4 | ||||||
|
|
||||||
| # Configuring case dictionary | ||||||
| print( | ||||||
|
There was a problem hiding this comment. Suggestion: Module-level side effect: the code prints the JSON at import time which causes unexpected output whenever this module is imported; wrap the printing in an Severity Level: Minor
Suggested change
Why it matters? ⭐The module currently performs a print at import time which is a side effect and can surprise callers who import this module. Wrapping the print in if name == "main": is the conventional, low-risk fix to ensure the JSON is only emitted when the script is executed directly. Prompt for AI Agent 🤖This is a comment left during a code review.
**Path:** examples/2D_mibm_cylinder_in_cross_flow/case.py
**Line:** 8:8
**Comment:**
*Logic Error: Module-level side effect: the code prints the JSON at import time which causes unexpected output whenever this module is imported; wrap the printing in an `if __name__ == "__main__":` guard so the JSON is only printed when the script is executed directly.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise. |
||||||
| json.dumps( | ||||||
| { | ||||||
| # Logistics | ||||||
| "run_time_info": "T", | ||||||
| # Computational Domain Parameters | ||||||
| # For these computations, the cylinder is placed at the (0,0,0) | ||||||
| # domain origin. | ||||||
| # axial direction | ||||||
| "x_domain%beg": 0.0e00, | ||||||
| "x_domain%end": 6.0e-03, | ||||||
| # r direction | ||||||
| "y_domain%beg": 0.0e00, | ||||||
| "y_domain%end": 6.0e-03, | ||||||
| "cyl_coord": "F", | ||||||
| "m": 200, | ||||||
| "n": 200, | ||||||
| "p": 0, | ||||||
| "dt": 6.0e-6, | ||||||
| "t_step_start": 0, | ||||||
| "t_step_stop": 10000, # 10000, | ||||||
| "t_step_save": 100, | ||||||
| # Simulation Algorithm Parameters | ||||||
| # Only one patches are necessary, the air tube | ||||||
| "num_patches": 1, | ||||||
| # Use the 5 equation model | ||||||
| "model_eqns": 2, | ||||||
| "alt_soundspeed": "F", | ||||||
| # One fluids: air | ||||||
| "num_fluids": 1, | ||||||
| # time step | ||||||
| "mpp_lim": "F", | ||||||
| # Correct errors when computing speed of sound | ||||||
| "mixture_err": "T", | ||||||
| # Use TVD RK3 for time marching | ||||||
| "time_stepper": 3, | ||||||
| # Use WENO5 | ||||||
| "weno_order": 5, | ||||||
| "weno_eps": 1.0e-16, | ||||||
| "weno_Re_flux": "T", | ||||||
| "weno_avg": "T", | ||||||
| "avg_state": 2, | ||||||
| "mapped_weno": "T", | ||||||
| "null_weights": "F", | ||||||
| "mp_weno": "T", | ||||||
| "riemann_solver": 2, | ||||||
| "wave_speeds": 1, | ||||||
| # We use ghost-cell | ||||||
| "bc_x%beg": -3, | ||||||
| "bc_x%end": -3, | ||||||
| "bc_y%beg": -3, | ||||||
| "bc_y%end": -3, | ||||||
| # Set IB to True and add 1 patch | ||||||
| "ib": "T", | ||||||
| "num_ibs": 1, | ||||||
| "viscous": "T", | ||||||
| # Formatted Database Files Structure Parameters | ||||||
| "format": 1, | ||||||
| "precision": 2, | ||||||
| "prim_vars_wrt": "T", | ||||||
| "E_wrt": "T", | ||||||
| "parallel_io": "T", | ||||||
| # Patch: Constant Tube filled with air | ||||||
| # Specify the cylindrical air tube grid geometry | ||||||
| "patch_icpp(1)%geometry": 3, | ||||||
| "patch_icpp(1)%x_centroid": 3.0e-03, | ||||||
| # Uniform medium density, centroid is at the center of the domain | ||||||
| "patch_icpp(1)%y_centroid": 3.0e-03, | ||||||
| "patch_icpp(1)%length_x": 6.0e-03, | ||||||
| "patch_icpp(1)%length_y": 6.0e-03, | ||||||
| # Specify the patch primitive variables | ||||||
| "patch_icpp(1)%vel(1)": 0.05e00, | ||||||
| "patch_icpp(1)%vel(2)": 0.0e00, | ||||||
| "patch_icpp(1)%pres": 1.0e00, | ||||||
| "patch_icpp(1)%alpha_rho(1)": 1.0e00, | ||||||
| "patch_icpp(1)%alpha(1)": 1.0e00, | ||||||
| # Patch: Cylinder Immersed Boundary | ||||||
| "patch_ib(1)%geometry": 2, | ||||||
| "patch_ib(1)%x_centroid": 1.5e-03, | ||||||
| "patch_ib(1)%y_centroid": 4.5e-03, | ||||||
| "patch_ib(1)%radius": 0.3e-03, | ||||||
| "patch_ib(1)%slip": "F", | ||||||
| "patch_ib(1)%moving_ibm": 2, | ||||||
| "patch_ib(1)%vel(2)": -0.1, | ||||||
| "patch_ib(1)%angles(1)": 0.0, # x-axis rotation in radians | ||||||
| "patch_ib(1)%angles(2)": 0.0, # y-axis rotation | ||||||
| "patch_ib(1)%angles(3)": 0.0, # z-axis rotation | ||||||
| "patch_ib(1)%angular_vel(1)": 0.0, # x-axis rotational velocity in radians per second | ||||||
| "patch_ib(1)%angular_vel(2)": 0.0, # y-axis rotation | ||||||
| "patch_ib(1)%angular_vel(3)": 0.0, # z-axis rotation | ||||||
| "patch_ib(1)%mass": 1.0e-6, # z-axis rotation | ||||||
|
Contributor
There was a problem hiding this comment. P2: Misleading comment: The comment says Prompt for AI agents
Suggested change
|
||||||
| # Fluids Physical Parameters | ||||||
| "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00), # 2.50(Not 1.40) | ||||||
| "fluid_pp(1)%pi_inf": 0, | ||||||
| "fluid_pp(1)%Re(1)": 2500000, | ||||||
| } | ||||||
| ) | ||||||
| ) | ||||||
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
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
| Original file line number | Diff line number | Diff line change | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|
|
|
@@ -472,6 +472,32 @@ contains | |||||||||
|
|
||||||||||
| end subroutine s_mpi_allreduce_sum | ||||||||||
|
|
||||||||||
| !> This subroutine follows the behavior of the s_mpi_allreduce_sum subroutine | ||||||||||
| !> with the additional feature that it reduces an array of vectors. | ||||||||||
| impure subroutine s_mpi_allreduce_vectors_sum(var_loc, var_glb, num_vectors, vector_length) | ||||||||||
|
|
||||||||||
| integer, intent(in) :: num_vectors, vector_length | ||||||||||
| real(wp), dimension(:, :), intent(in) :: var_loc | ||||||||||
| real(wp), dimension(:, :), intent(out) :: var_glb | ||||||||||
|
Comment on lines
+479
to
+481
There was a problem hiding this comment. Suggestion: Intent mismatch for in-place reduction: Severity Level: Minor
Suggested change
Why it matters? ⭐This is a valid correctness issue. INTENT(OUT) permits the compiler to treat the actual argument as undefined on entry (and even reallocate/erase it), which makes using MPI_IN_PLACE unsafe when the caller passes the same buffer for send and receive. Changing var_glb to INTENT(INOUT) accurately documents and guarantees the value is preserved on entry and is the correct contract for potential in-place MPI calls. Prompt for AI Agent 🤖This is a comment left during a code review.
**Path:** src/common/m_mpi_common.fpp
**Line:** 479:481
**Comment:**
*Logic Error: Intent mismatch for in-place reduction: `var_glb` is declared with INTENT(OUT) but the subroutine uses MPI_IN_PLACE when the caller passes the same buffer for send and receive; INTENT(OUT) allows the compiler to undefine or reallocate the actual argument before the call which makes in-place MPI operations unsafe. Change `var_glb` to INTENT(INOUT) so the caller's data is preserved on entry for MPI_IN_PLACE use.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise. |
||||||||||
|
|
||||||||||
| #ifdef MFC_MPI | ||||||||||
| integer :: ierr !< Generic flag used to identify and report MPI errors | ||||||||||
|
|
||||||||||
| ! Performing the reduction procedure | ||||||||||
| if (loc(var_loc) == loc(var_glb)) then | ||||||||||
| call MPI_Allreduce(MPI_IN_PLACE, var_glb, num_vectors*vector_length, & | ||||||||||
| mpi_p, MPI_SUM, MPI_COMM_WORLD, ierr) | ||||||||||
| else | ||||||||||
| call MPI_Allreduce(var_loc, var_glb, num_vectors*vector_length, & | ||||||||||
| mpi_p, MPI_SUM, MPI_COMM_WORLD, ierr) | ||||||||||
| end if | ||||||||||
|
|
||||||||||
| #else | ||||||||||
| var_glb(1:num_vectors, 1:vector_length) = var_loc(1:num_vectors, 1:vector_length) | ||||||||||
| #endif | ||||||||||
|
Contributor
There was a problem hiding this comment. P2: Missing non-MPI fallback: when Prompt for AI agents
Suggested change
✅ Addressed in |
||||||||||
|
|
||||||||||
| end subroutine s_mpi_allreduce_vectors_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 | ||||||||||
|
|
||||||||||
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
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -560,6 +560,8 @@ contains | |
| patch_ib(i)%vel(:) = 0._wp | ||
| patch_ib(i)%angles(:) = 0._wp | ||
| patch_ib(i)%angular_vel(:) = 0._wp | ||
| patch_ib(i)%mass = dflt_real | ||
| patch_ib(i)%moment = dflt_real | ||
|
|
||
| ! sets values of a rotation matrix which can be used when calculating rotations | ||
| patch_ib(i)%rotation_matrix = 0._wp | ||
|
|
||
|
Contributor
There was a problem hiding this comment. High-level SuggestionThe force and torque calculations in Solution Walkthrough:Before:subroutine s_compute_ib_forces(pressure)
...
do i = 0, m
do j = 0, n
do k = 0, p
...
! Incorrect pressure gradient using cell-center coordinates in denominator
pressure_divergence(1) = (pressure(i+1,j,k) - pressure(i-1,j,k)) / (2._wp * x_cc(i))
pressure_divergence(2) = (pressure(i,j+1,k) - pressure(i,j-1,k)) / (2._wp * y_cc(j))
! Incorrect cell volume calculation using coordinates
cell_volume = x_cc(i) * y_cc(j)
if (num_dims == 3) then
pressure_divergence(3) = (pressure(i,j,k+1) - pressure(i,j,k-1)) / (2._wp * z_cc(k))
cell_volume = cell_volume * z_cc(k)
end if
forces(ib_idx, :) = forces(ib_idx, :) - (pressure_divergence * cell_volume)
...
end do
end do
end do
end subroutine
After:subroutine s_compute_ib_forces(pressure)
...
! Grid spacings (dx, dy, dz) should be available or computed
! e.g., dx = x_cc(1) - x_cc(0) for a uniform grid
do i = 0, m
do j = 0, n
do k = 0, p
...
! Correct pressure gradient using grid spacing
pressure_divergence(1) = (pressure(i+1,j,k) - pressure(i-1,j,k)) / (2._wp * dx)
pressure_divergence(2) = (pressure(i,j+1,k) - pressure(i,j-1,k)) / (2._wp * dy)
! Correct cell volume calculation
cell_volume = dx * dy
if (num_dims == 3) then
pressure_divergence(3) = (pressure(i,j,k+1) - pressure(i,j,k-1)) / (2._wp * dz)
cell_volume = cell_volume * dz
end if
forces(ib_idx, :) = forces(ib_idx, :) - (pressure_divergence * cell_volume)
...
end do
end do
end do
end subroutine
|
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
| Original file line number | Diff line number | Diff line change | ||
|---|---|---|---|---|
|
|
@@ -101,7 +101,6 @@ contains | |||
| end if | ||||
| call s_update_ib_rotation_matrix(i) | ||||
| end do | ||||
| $:GPU_ENTER_DATA(copyin='[patch_ib]') | ||||
|
|
||||
| ! Allocating the patch identities bookkeeping variable | ||||
|
|
@@ -197,6 +196,19 @@ contains | |||
| type(ghost_point) :: gp | ||||
| type(ghost_point) :: innerp | ||||
|
|
||||
| ! set the Moving IBM interior Pressure Values | ||||
| $:GPU_PARALLEL_LOOP(private='[i,j,k]', copyin='[E_idx]', collapse=3) | ||||
| do l = 0, p | ||||
| do k = 0, n | ||||
| do j = 0, m | ||||
| if (ib_markers%sf(j, k, l) /= 0) then | ||||
| q_prim_vf(E_idx)%sf(j, k, l) = 1._wp | ||||
| end if | ||||
| end do | ||||
| end do | ||||
| end do | ||||
| $:END_GPU_PARALLEL_LOOP() | ||||
|
|
||||
| if (num_gps > 0) then | ||||
| $:GPU_PARALLEL_LOOP(private='[i,physical_loc,dyn_pres,alpha_rho_IP, alpha_IP,pres_IP,vel_IP,vel_g,vel_norm_IP,r_IP, v_IP,pb_IP,mv_IP,nmom_IP,presb_IP,massv_IP,rho, gamma,pi_inf,Re_K,G_K,Gs,gp,innerp,norm,buf, radial_vector, rotation_velocity, j,k,l,q,qv_K,c_IP,nbub,patch_id]') | ||||
| do i = 1, num_gps | ||||
|
|
@@ -244,6 +256,18 @@ contains | |||
| if (surface_tension) then | ||||
| q_prim_vf(c_idx)%sf(j, k, l) = c_IP | ||||
| end if | ||||
|
|
||||
| ! set the pressure | ||||
| if (patch_ib(patch_id)%moving_ibm == 0) then | ||||
| q_prim_vf(E_idx)%sf(j, k, l) = pres_IP | ||||
| else | ||||
| q_prim_vf(E_idx)%sf(j, k, l) = 0._wp | ||||
| $:GPU_LOOP(parallelism='[seq]') | ||||
| do q = 1, num_fluids | ||||
| q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) + pres_IP/(1._wp - 2._wp*abs(levelset%sf(j, k, l, patch_id)*alpha_rho_IP(q)/pres_IP)*dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, levelset_norm%sf(j, k, l, patch_id, :))) | ||||
| end do | ||||
| end if | ||||
|
|
||||
| if (model_eqns /= 4) then | ||||
| ! If in simulation, use acc mixture subroutines | ||||
| if (elasticity) then | ||||
|
|
@@ -969,6 +993,78 @@ contains | |||
|
|
||||
| end subroutine s_update_mib | ||||
|
|
||||
| ! compute the surface integrals of the IB via a volume integraion method described in | ||||
| ! "A coupled IBM/Euler-Lagrange framework for simulating shock-induced particle size segregation" | ||||
| ! by Archana Sridhar and Jesse Capecelatro | ||||
| subroutine s_compute_ib_forces(pressure) | ||||
|
|
||||
| real(wp), dimension(0:m, 0:n, 0:p), intent(in) :: pressure | ||||
|
|
||||
| integer :: i, j, k, ib_idx | ||||
| real(wp), dimension(num_ibs, 3) :: forces, torques | ||||
| real(wp), dimension(1:3) :: pressure_divergence, radial_vector | ||||
| real(wp) :: cell_volume, dx, dy, dz | ||||
|
|
||||
| forces = 0._wp | ||||
| torques = 0._wp | ||||
|
|
||||
| ! TODO :: This is currently only valid inviscid, and needs to be extended to add viscocity | ||||
| $:GPU_PARALLEL_LOOP(private='[ib_idx,radial_vector,pressure_divergence,cell_volume, dx, dy, dz]', copy='[forces,torques]', copyin='[ib_markers]', collapse=3) | ||||
| do i = 0, m | ||||
| do j = 0, n | ||||
| do k = 0, p | ||||
| ib_idx = ib_markers%sf(i, j, k) | ||||
| if (ib_idx /= 0) then ! only need to compute the gradient for cells inside a IB | ||||
| if (patch_ib(ib_idx)%moving_ibm == 2) then ! make sure that this IB has 2-way coupling enabled | ||||
| ! get the vector pointing to the grid cell from the IB centroid | ||||
| if (num_dims == 3) then | ||||
| radial_vector = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid] | ||||
| else | ||||
| radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, patch_ib(ib_idx)%y_centroid, 0._wp] | ||||
| end if | ||||
| dx = x_cc(i + 1) - x_cc(i) | ||||
| dy = y_cc(j + 1) - y_cc(j) | ||||
|
|
||||
| ! use a finite difference to compute the 2D components of the gradient of the pressure and cell volume | ||||
| pressure_divergence(1) = (pressure(i + 1, j, k) - pressure(i - 1, j, k))/(2._wp*dx) | ||||
| pressure_divergence(2) = (pressure(i, j + 1, k) - pressure(i, j - 1, k))/(2._wp*dy) | ||||
| cell_volume = dx*dy | ||||
|
|
||||
| ! add the 3D component, if we are working in 3 dimensions | ||||
| if (num_dims == 3) then | ||||
| dz = z_cc(k + 1) - z_cc(k) | ||||
| pressure_divergence(3) = (pressure(i, j, k + 1) - pressure(i, j, k - 1))/(2._wp*dz) | ||||
| cell_volume = cell_volume*dz | ||||
| else | ||||
| pressure_divergence(3) = 0._wp | ||||
| end if | ||||
|
|
||||
| ! Update the force values atomically to prevent race conditions | ||||
| $:GPU_ATOMIC(atomic='update') | ||||
| forces(ib_idx, :) = forces(ib_idx, :) - (pressure_divergence*cell_volume) | ||||
| $:GPU_ATOMIC(atomic='update') | ||||
| torques(ib_idx, :) = torques(ib_idx, :) - (cross_product(radial_vector, pressure_divergence)*cell_volume) | ||||
| end if | ||||
| end if | ||||
| end do | ||||
| end do | ||||
| end do | ||||
| $:END_GPU_PARALLEL_LOOP() | ||||
|
|
||||
| ! reduce the forces across all MPI ranks | ||||
| call s_mpi_allreduce_vectors_sum(forces, forces, num_ibs, 3) | ||||
| call s_mpi_allreduce_vectors_sum(torques, torques, num_ibs, 3) | ||||
|
|
||||
| ! apply the summed forces | ||||
| do i = 1, num_ibs | ||||
| patch_ib(i)%force(:) = forces(i, :) | ||||
| patch_ib(i)%torque(:) = matmul(patch_ib(i)%rotation_matrix_inverse, torques(i, :)) ! torques must be computed in the local coordinates of the IB | ||||
| end do | ||||
|
|
||||
| print *, forces(1, 1:2) | ||||
|
Contributor
There was a problem hiding this comment. Remove debug print statement. This debug output will produce excessive noise in production runs. Remove or guard with a debug flag. - print *, forces(1, 1:2)📝 Committable suggestion
Suggested change
🤖 Prompt for AI Agents
Contributor
There was a problem hiding this comment. P2: Debug print statement should be removed before merging. This will pollute output logs and may cause performance issues in parallel execution. Prompt for AI agents |
||||
|
|
||||
| end subroutine s_compute_ib_forces | ||||
|
|
||||
| !> Subroutine to deallocate memory reserved for the IBM module | ||||
| impure subroutine s_finalize_ibm_module() | ||||
|
|
||||
|
|
@@ -978,6 +1074,77 @@ contains | |||
|
|
||||
| end subroutine s_finalize_ibm_module | ||||
|
|
||||
| subroutine s_compute_moment_of_inertia(ib_marker, axis) | ||||
|
|
||||
| real(wp), dimension(3), optional :: axis !< the axis about which we compute the moment. Only required in 3D. | ||||
| integer, intent(in) :: ib_marker | ||||
|
|
||||
| real(wp) :: moment, distance_to_axis, cell_volume | ||||
| real(wp), dimension(3) :: position, closest_point_along_axis, vector_to_axis | ||||
| integer :: i, j, k, count | ||||
|
|
||||
| if (p == 0) then | ||||
| axis = [0, 1, 0] | ||||
| else if (sqrt(sum(axis**2)) == 0) then | ||||
| ! if the object is not actually rotating at this time, return a dummy value and exit | ||||
| patch_ib(ib_marker)%moment = 1._wp | ||||
| return | ||||
| else | ||||
| axis = axis/sqrt(sum(axis)) | ||||
| end if | ||||
|
Comment on lines
+1087
to
+1095
Contributor
There was a problem hiding this comment. Optional argument handling issue in 2D case. When -subroutine s_compute_moment_of_inertia(ib_marker, axis)
+subroutine s_compute_moment_of_inertia(ib_marker, axis_in)
- real(wp), dimension(3), optional :: axis !< the axis about which we compute the moment. Only required in 3D.
+ real(wp), dimension(3), intent(in), optional :: axis_in !< the axis about which we compute the moment. Only required in 3D.
integer, intent(in) :: ib_marker
real(wp) :: moment, distance_to_axis, cell_volume
- real(wp), dimension(3) :: position, closest_point_along_axis, vector_to_axis
+ real(wp), dimension(3) :: position, closest_point_along_axis, vector_to_axis, axis
integer :: i, j, k, count
if (p == 0) then
- axis = [0, 1, 0]
- else if (sqrt(sum(axis**2)) == 0) then
+ axis = [0._wp, 0._wp, 1._wp]
+ else if (.not. present(axis_in) .or. sqrt(sum(axis_in**2)) < 1.e-16_wp) then
! if the object is not actually rotating at this time, return a dummy value and exit
patch_ib(ib_marker)%moment = 1._wp
return
else
- axis = axis/sqrt(sum(axis))
+ axis = axis_in/sqrt(sum(axis_in**2))
end if
🤖 Prompt for AI Agents |
||||
|
|
||||
| ! if the IB is in 2D or a 3D sphere, we can compute this exactly | ||||
| if (patch_ib(ib_marker)%geometry == 2) then ! circle | ||||
| patch_ib(ib_marker)%moment = 0.5*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%radius)**2 | ||||
| elseif (patch_ib(ib_marker)%geometry == 3) then ! rectangle | ||||
| patch_ib(ib_marker)%moment = patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%length_x**2 + patch_ib(ib_marker)%length_y**2)/6._wp | ||||
| elseif (patch_ib(ib_marker)%geometry == 8) then ! sphere | ||||
| patch_ib(ib_marker)%moment = 0.4*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%radius)**2 | ||||
|
|
||||
| else ! we do not have an analytic moment of inertia calculation and need to approximate it directly | ||||
| count = 0 | ||||
| moment = 0._wp | ||||
| cell_volume = (x_cc(1) - x_cc(0))*(y_cc(1) - y_cc(0)) ! computed without grid stretching. Update in the loop to perform with stretching | ||||
| if (p /= 0) then | ||||
| cell_volume = cell_volume*(z_cc(1) - z_cc(0)) | ||||
| end if | ||||
|
|
||||
| $:GPU_PARALLEL_LOOP(private='[position,closest_point_along_axis,vector_to_axis,distance_to_axis]', copy='[moment,count]', copyin='[ib_marker,cell_volume,axis]', collapse=3) | ||||
| do i = 0, m | ||||
| do j = 0, n | ||||
| do k = 0, p | ||||
| if (ib_markers%sf(i, j, k) == ib_marker) then | ||||
| $:GPU_ATOMIC(atomic='update') | ||||
| count = count + 1 ! increment the count of total cells in the boundary | ||||
|
|
||||
| ! get the position in local coordinates so that the axis passes through 0, 0, 0 | ||||
| if (p == 0) then | ||||
| position = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, 0._wp] | ||||
| else | ||||
| position = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%z_centroid] | ||||
| end if | ||||
|
|
||||
| ! project the position along the axis to find the closest distance to the rotation axis | ||||
| closest_point_along_axis = axis*dot_product(axis, position) | ||||
| vector_to_axis = position - closest_point_along_axis | ||||
| distance_to_axis = dot_product(vector_to_axis, vector_to_axis) ! saves the distance to the axis squared | ||||
|
|
||||
| ! compute the position component of the moment | ||||
| $:GPU_ATOMIC(atomic='update') | ||||
| moment = moment + distance_to_axis | ||||
| end if | ||||
| end do | ||||
| end do | ||||
| end do | ||||
| $:END_GPU_PARALLEL_LOOP() | ||||
|
|
||||
| ! write the final moment assuming the points are all uniform density | ||||
| patch_ib(ib_marker)%moment = moment*patch_ib(ib_marker)%mass/(count*cell_volume) | ||||
| $:GPU_UPDATE(device='[patch_ib(ib_marker)%moment]') | ||||
| end if | ||||
|
|
||||
| end subroutine s_compute_moment_of_inertia | ||||
|
|
||||
| function cross_product(a, b) result(c) | ||||
| implicit none | ||||
| real(wp), intent(in) :: a(3), b(3) | ||||
|
|
||||
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
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -203,7 +203,7 @@ contains | |
|
|
||
| do i = 1, num_ibs | ||
| #:for VAR in [ 'radius', 'length_x', 'length_y', 'length_z', & | ||
| & 'x_centroid', 'y_centroid', 'z_centroid', 'c', 'm', 'p', 't', 'theta', 'slip', 'mass'] | ||
| call MPI_BCAST(patch_ib(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) | ||
| #:endfor | ||
| #:for VAR in ['vel', 'angular_vel', 'angles'] | ||
|
|
||
Oops, something went wrong.
Oops, something went wrong.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Suggestion: Unused import:
mathis imported but never used anywhere; keep dead imports out of the codebase to avoid confusion and maintain clean dependencies — remove the import. [possible bug]Severity Level: Critical 🚨
Why it matters? ⭐
The module imports math but never uses it anywhere in the current PR diff. Removing the import removes dead code and linter warnings and has no functional impact.
Prompt for AI Agent 🤖