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

Draft
SVS87 wants to merge 8 commits into
base: master
Choose a base branch
from
Draft

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

Show file tree
Hide file tree
Changes from all 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_fluid_rectangle/case.py
View file
Open in desktop

@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.

Not sure about this name, as it could get confused with IB patch code. Maybe speciify it is a fluid patch?

@SVS87 SVS87 Jun 21, 2026

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

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

Changing it to 2D_rotated_fluid_rectangle. Would that work?

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,
}
)
)
5 changes: 5 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 Expand Up @@ -333,6 +337,7 @@ module m_derived_types
real(wp), dimension(1:3) :: model_translate !< Translation of the STL object.
real(wp), dimension(1:3) :: model_scale !< Scale factor for the STL object.
real(wp) :: model_threshold !< Threshold to turn on smooth STL patch.
real(wp), dimension(1:3) :: model_rotate !< Rotation of the STL object (radians, per axis).
end type ib_stl_parameters

type ib_patch_parameters
Expand Down
2 changes: 1 addition & 1 deletion src/common/m_model.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -887,7 +887,7 @@ contains
model = f_model_read(stl_models(stl_id)%model_filepath)
params%scale(:) = stl_models(stl_id)%model_scale(:)
params%translate(:) = stl_models(stl_id)%model_translate(:)
params%rotate(:) = 0._wp
params%rotate(:) = stl_models(stl_id)%model_rotate(:)
params%spc = num_ray
params%threshold = stl_models(stl_id)%model_threshold

Expand Down
47 changes: 46 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,48 @@ 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


! construct the x, y, and z rotation matrices

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

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

! 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
rotation(3, 1,:) = [cos(angles(3)), -sin(angles(3)), 0._wp]
rotation(3, 2,:) = [sin(angles(3)), cos(angles(3)), 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
2 changes: 2 additions & 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 Expand Up @@ -336,6 +337,7 @@ contains
stl_models(i)%model_filepath(:) = dflt_char
stl_models(i)%model_translate(:) = 0._wp
stl_models(i)%model_scale(:) = 1._wp
stl_models(i)%model_rotate(:) = 0._wp
stl_models(i)%model_threshold = ray_tracing_threshold
end do

Expand Down
Loading
Loading