Skip to content

feature: integrate rotations into ICPP patches (#1594) #1612

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Closed
SVS87 wants to merge 8 commits into from
Closed

feature: integrate rotations into ICPP patches (#1594) #1612

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
c4dea86
feature: integrate rotations into ICPP patches (#1594)
SVS87 Jun 20, 2026
196a73f
Merge branch 'master' into feature/icpp-patch-rotations
SVS87 Jun 20, 2026
f74e83c
Address review: extend ICPP rotation to all geometries, remove IB rot…
SVS87 Jun 21, 2026
7372c5d
Merge branch 'feature/icpp-patch-rotations' of https://github.com/SVS…
SVS87 Jun 21, 2026
ae7c8b5
Add model_rotate for STL patches; complete rotation rollout
SVS87 Jun 21, 2026
db32925
Add golden files for 2D_rotated_fluid_rectangle example
SVS87 Jun 21, 2026
2b2d4ac
Re-trigger CI
SVS87 Jun 21, 2026
cf6be2f
Merge branch 'master' into feature/icpp-patch-rotations
SVS87 Jun 22, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 84 additions & 0 deletions examples/2D_rotated_rectangle/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
#!/usr/bin/env python3
import json

# Configuring case dictionary
print(
json.dumps(
{
# Logistics
"run_time_info": "T",
# Computational Domain Parameters
"x_domain%beg": 0.0e00,
"x_domain%end": 1.0e00,
"y_domain%beg": 0.0e00,
"y_domain%end": 1.0e00,
"m": 99,
"n": 99,
"p": 0,
"dt": 1.0e-06,
"t_step_start": 0,
"t_step_stop": 50,
"t_step_save": 10,
# Simulation Algorithm Parameters
"num_patches": 2,
"model_eqns": "6eq",
"alt_soundspeed": "F",
"num_fluids": 2,
"mpp_lim": "T",
"mixture_err": "T",
"time_stepper": "rk3",
"weno_order": 5,
"weno_eps": 1.0e-16,
"weno_Re_flux": "F",
"weno_avg": "F",
"mapped_weno": "T",
"null_weights": "F",
"mp_weno": "F",
"riemann_solver": "hllc",
"wave_speeds": "direct",
"avg_state": "arithmetic",
"bc_x%beg": -3,
"bc_x%end": -3,
"bc_y%beg": -3,
"bc_y%end": -3,
# Formatted Database Files Structure Parameters
"format": "silo",
"precision": "double",
"prim_vars_wrt": "T",
"parallel_io": "T",
# Patch 1: Base
"patch_icpp(1)%geometry": 3,
"patch_icpp(1)%x_centroid": 0.5e00,
"patch_icpp(1)%y_centroid": 0.5e00,
"patch_icpp(1)%length_x": 1.0e00,
"patch_icpp(1)%length_y": 1.0e00,
"patch_icpp(1)%vel(1)": 0.0e00,
"patch_icpp(1)%vel(2)": 0.0e00,
"patch_icpp(1)%pres": 1.0e05,
"patch_icpp(1)%alpha_rho(1)": 1.0e00,
"patch_icpp(1)%alpha_rho(2)": 1.0e-12,
"patch_icpp(1)%alpha(1)": 1.0 - 1.0e-12,
"patch_icpp(1)%alpha(2)": 1.0e-12,
# Patch 2: Rotated rectangle (45 degrees)
"patch_icpp(2)%geometry": 3,
"patch_icpp(2)%x_centroid": 0.5e00,
"patch_icpp(2)%y_centroid": 0.5e00,
"patch_icpp(2)%length_x": 0.4e00,
"patch_icpp(2)%length_y": 0.15e00,
"patch_icpp(2)%angles(3)": 0.785398e00,
"patch_icpp(2)%alter_patch(1)": "T",
"patch_icpp(2)%vel(1)": 0.0e00,
"patch_icpp(2)%vel(2)": 0.0e00,
"patch_icpp(2)%pres": 1.0e05,
"patch_icpp(2)%alpha_rho(1)": 1.0e-12,
"patch_icpp(2)%alpha_rho(2)": 1.0e00,
"patch_icpp(2)%alpha(1)": 1.0e-12,
"patch_icpp(2)%alpha(2)": 1.0 - 1.0e-12,
# Fluids Physical Parameters
"fluid_pp(1)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
"fluid_pp(1)%pi_inf": 0.0e00,
"fluid_pp(2)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
"fluid_pp(2)%pi_inf": 0.0e00,
}
)
)
4 changes: 4 additions & 0 deletions src/common/m_derived_types.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -284,6 +284,10 @@ module m_derived_types
! Geometry 14 (3D spherical harmonic): sph_har_coeff(l,m) for real Y_lm
real(wp), dimension(0:max_sph_harm_degree,-max_sph_harm_degree:max_sph_harm_degree) :: sph_har_coeff
real(wp), dimension(3) :: normal !< Patch orientation normal vector (x, y, z)
! Patch rotation (shared framework; see m_patch_geometries)
real(wp), dimension(1:3) :: angles !< Rotation angles about x, y, z axes
real(wp), dimension(1:3,1:3) :: rotation_matrix !< Rotation matrix built from angles
real(wp), dimension(1:3,1:3) :: rotation_matrix_inverse !< Inverse of rotation matrix
logical, dimension(0:num_patches_max - 1) :: alter_patch !< Overwrite permissions for preceding patches
logical :: smoothen !< Whether patch boundaries are smoothed across cells
integer :: smooth_patch_id !< Identity (id) of the patch with which current patch is to get smoothed
Expand Down
51 changes: 50 additions & 1 deletion src/common/m_patch_geometries.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,8 @@ module m_patch_geometries

implicit none

public :: f_is_inside_sphere, f_is_inside_cylinder, f_is_inside_cuboid, f_is_inside_airfoil, f_is_inside_ellipse
public :: f_is_inside_sphere, f_is_inside_cylinder, f_is_inside_cuboid, f_is_inside_airfoil, f_is_inside_ellipse, &
& s_compute_rotation_matrix

contains

Expand Down Expand Up @@ -140,4 +141,52 @@ contains

end function f_is_inside_ellipse

!> Compute a rotation matrix for converting to the rotating frame of the boundary
subroutine s_compute_rotation_matrix(angles, rotation_matrix, rotation_matrix_inverse)

$:GPU_ROUTINE(parallelism='[seq]')

real(wp), dimension(1:3), intent(in) :: angles
real(wp), dimension(1:3,1:3), intent(out) :: rotation_matrix
real(wp), dimension(1:3,1:3), intent(out) :: rotation_matrix_inverse
real(wp), dimension(3, 3, 3) :: rotation

@danieljvickers danieljvickers Jun 20, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Making this 3x3x3 is a nice touch

real(wp) :: angle

@danieljvickers danieljvickers Jun 20, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Is there a reason to add this additional angle value?

@SVS87 SVS87 Jun 21, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Removed it, it was unnecessary. Updated change just uses angles(1, 2, 3) respectively


! construct the x, y, and z rotation matrices

if (num_dims == 3) then
! also compute the x and y axes in 3D
angle = angles(1)
rotation(1, 1,:) = [1._wp, 0._wp, 0._wp]
rotation(1, 2,:) = [0._wp, cos(angle), -sin(angle)]
rotation(1, 3,:) = [0._wp, sin(angle), cos(angle)]

angle = angles(2)
rotation(2, 1,:) = [cos(angle), 0._wp, sin(angle)]
rotation(2, 2,:) = [0._wp, 1._wp, 0._wp]
rotation(2, 3,:) = [-sin(angle), 0._wp, cos(angle)]

! apply the y rotation to the x rotation
rotation_matrix(:,:) = matmul(rotation(1,:,:), rotation(2,:,:))
rotation_matrix_inverse(:,:) = matmul(transpose(rotation(2,:,:)), transpose(rotation(1,:,:)))
end if

! z component first, since it applies in 2D and 3D
angle = angles(3)
rotation(3, 1,:) = [cos(angle), -sin(angle), 0._wp]
rotation(3, 2,:) = [sin(angle), cos(angle), 0._wp]
rotation(3, 3,:) = [0._wp, 0._wp, 1._wp]

if (num_dims == 3) then
! apply the z rotation to the xy rotation in 3D
rotation_matrix(:,:) = matmul(rotation_matrix(:,:), rotation(3,:,:))
rotation_matrix_inverse(:,:) = matmul(transpose(rotation(3,:,:)), rotation_matrix_inverse(:,:))
else
! write out only the z rotation in 2D
rotation_matrix(:,:) = rotation(3,:,:)
rotation_matrix_inverse(:,:) = transpose(rotation(3,:,:))
end if

end subroutine s_compute_rotation_matrix

end module m_patch_geometries
1 change: 1 addition & 0 deletions src/pre_process/m_global_parameters.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,6 +205,7 @@ contains
patch_icpp(i)%epsilon = dflt_real
patch_icpp(i)%beta = dflt_real
patch_icpp(i)%normal = dflt_real
patch_icpp(i)%angles = 0._wp

@sbryngelson sbryngelson Jun 23, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

This 0 default is exactly the value non-root MPI ranks will keep, because patch_icpp%angles is never broadcast. In src/pre_process/m_mpi_proxy.fpp the patch_icpp array-field broadcast loop (around line 103) is #:for VAR in ['normal','radii','vel','tau_e','alpha_rho','alpha','fourier_cos','fourier_sin'] — angles is missing. patch_icpp is intentionally part of the manual broadcast residue (see toolchain/.../fortran_gen.py, which lists patch_icpp as manual), so it is NOT covered by generated_bcast.fpp; the manual loop is the only broadcaster.

Consequence: only rank 0 reads angles from the case file. Each rank then computes rotation_matrix(_inverse) locally from its own angles, so non-root ranks build the identity → their subdomains get the patch unrotated. In a multi-rank pre_process run this produces a torn initial condition (rotated on rank 0, axis-aligned elsewhere), silently and dependent on domain decomposition. Affects every rotated geometry in this PR.

Fix: add 'angles' to that #:for VAR list in m_mpi_proxy.fpp (rotation_matrix/rotation_matrix_inverse are recomputed per-rank and don't need broadcasting). Also regenerate the golden files afterward, since any that were produced multi-rank would encode the pre-fix torn output.

@sbryngelson sbryngelson Jun 23, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Concrete fix (GitHub won't let me post this as a click-to-apply suggestion because the target line in src/pre_process/m_mpi_proxy.fpp — the patch_icpp array-field broadcast loop, ~line 103 — isn't part of this PR's diff). The one-line change is to add 'angles' to that loop:

            #:for VAR in [ 'angles', 'normal', 'radii', 'vel', 'tau_e', 'alpha_rho', 'alpha', &
                'fourier_cos', 'fourier_sin' ]
                call MPI_BCAST(patch_icpp(i)%${VAR}$, size(patch_icpp(i)%${VAR}$), mpi_p, 0, MPI_COMM_WORLD, ierr)
            #:endfor

size(patch_icpp(i)%angles) == 3, matching the other rank-1 members already in the loop, so no special-casing is needed. rotation_matrix/rotation_matrix_inverse are recomputed per-rank and don't need broadcasting. After applying, regenerate any goldens that were produced multi-rank.

patch_icpp(i)%radii = dflt_real
patch_icpp(i)%alter_patch = .false.
patch_icpp(i)%alter_patch(0) = .true.
Expand Down
16 changes: 14 additions & 2 deletions src/pre_process/m_icpp_patches.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,9 @@ contains
print *, 'Processing patch', i
end if

call s_compute_rotation_matrix(patch_icpp(i)%angles, patch_icpp(i)%rotation_matrix, &
& patch_icpp(i)%rotation_matrix_inverse)

!> ICPP Patches
!> @{
! Spherical patch
Expand Down Expand Up @@ -103,6 +106,9 @@ contains
print *, 'Processing patch', i
end if

call s_compute_rotation_matrix(patch_icpp(i)%angles, patch_icpp(i)%rotation_matrix, &
& patch_icpp(i)%rotation_matrix_inverse)

!> ICPP Patches
!> @{
! Circular patch
Expand Down Expand Up @@ -147,6 +153,9 @@ contains
print *, 'Processing patch', i
end if

call s_compute_rotation_matrix(patch_icpp(i)%angles, patch_icpp(i)%rotation_matrix, &
& patch_icpp(i)%rotation_matrix_inverse)

! Line segment patch
if (patch_icpp(i)%geometry == 1) then
call s_icpp_line_segment(i, patch_id_fp, q_prim_vf)
Expand Down Expand Up @@ -604,7 +613,7 @@ contains
type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf
integer :: i, j, k !< generic loop iterators
real(wp) :: pi_inf, gamma, lit_gamma !< Equation of state parameters

real(wp), dimension(1:3) :: xyz_local
@:HardcodedDimensionsExtrusion()
@:Hardcoded2DVariables()

Expand All @@ -630,7 +639,10 @@ contains
! Assign patch vars if cell is covered and patch has write permission
do j = 0, n
do i = 0, m
if (f_is_inside_cuboid(x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp, [length_x, length_y, 0._wp])) then

@danieljvickers danieljvickers Jun 20, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

I see where you added cuboids. I assume you are in the process of adding this to the remaining geometries? That is fine if so. It just needs to be integrated with all of them before merge.

! Rotate the centroid-relative grid point into the patch frame
xyz_local = matmul(patch_icpp(patch_id)%rotation_matrix_inverse, [x_cc(i) - x_centroid, y_cc(j) - y_centroid, &
& 0._wp])
if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), 0._wp, [length_x, length_y, 0._wp])) then
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)

Expand Down
41 changes: 3 additions & 38 deletions src/simulation/m_ib_patches.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -415,45 +415,10 @@ contains

$:GPU_ROUTINE(parallelism='[seq]')

integer, intent(in) :: patch_id
real(wp), dimension(3, 3, 3) :: rotation
real(wp) :: angle
integer, intent(in) :: patch_id

! construct the x, y, and z rotation matrices

if (num_dims == 3) then
! also compute the x and y axes in 3D
angle = patch_ib(patch_id)%angles(1)
rotation(1, 1,:) = [1._wp, 0._wp, 0._wp]
rotation(1, 2,:) = [0._wp, cos(angle), -sin(angle)]
rotation(1, 3,:) = [0._wp, sin(angle), cos(angle)]

angle = patch_ib(patch_id)%angles(2)
rotation(2, 1,:) = [cos(angle), 0._wp, sin(angle)]
rotation(2, 2,:) = [0._wp, 1._wp, 0._wp]
rotation(2, 3,:) = [-sin(angle), 0._wp, cos(angle)]

! apply the y rotation to the x rotation
patch_ib(patch_id)%rotation_matrix(:,:) = matmul(rotation(1,:,:), rotation(2,:,:))
patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(2,:,:)), transpose(rotation(1,:,:)))
end if

! z component first, since it applies in 2D and 3D
angle = patch_ib(patch_id)%angles(3)
rotation(3, 1,:) = [cos(angle), -sin(angle), 0._wp]
rotation(3, 2,:) = [sin(angle), cos(angle), 0._wp]
rotation(3, 3,:) = [0._wp, 0._wp, 1._wp]

if (num_dims == 3) then
! apply the z rotation to the xy rotation in 3D
patch_ib(patch_id)%rotation_matrix(:,:) = matmul(patch_ib(patch_id)%rotation_matrix(:,:), rotation(3,:,:))
patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(3,:,:)), &
& patch_ib(patch_id)%rotation_matrix_inverse(:,:))
else
! write out only the z rotation in 2D
patch_ib(patch_id)%rotation_matrix(:,:) = rotation(3,:,:)
patch_ib(patch_id)%rotation_matrix_inverse(:,:) = transpose(rotation(3,:,:))
end if
call s_compute_rotation_matrix(patch_ib(patch_id)%angles, patch_ib(patch_id)%rotation_matrix, &
& patch_ib(patch_id)%rotation_matrix_inverse)

end subroutine s_update_ib_rotation_matrix

Expand Down
Loading
Loading