feature: integrate rotations into ICPP patches (#1594)

Add a rotation framework to initial-condition (ICPP) patches, mirroring
the existing immersed-boundary rotation support.

- Extract the rotation-matrix computation into a new shared module
  m_patch_geometries, taking angle/matrix arrays as arguments so both
  ICPP and IB patches can call it (m_ib_patches now uses a thin wrapper).
- Add angles + rotation_matrix(_inverse) fields to ic_patch_parameters,
  defaulting angles to zero (identity rotation) so existing cases are
  unaffected.
- Apply the inverse-rotation transform in the rectangle inside-test as a
  proof of concept; remaining geometries are straightforward follow-up.
- Register patch_icpp(i)%angles in the parameter definitions so the angle
  is settable from a case file.
- Add a 2D_rotated_rectangle example demonstrating a 45-degree patch.

Author: SVS87

examples/2D_rotated_rectangle/case.py

@@ -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": 5.0e-07,
+            "t_step_start": 0,
+            "t_step_stop": 1000,
+            "t_step_save": 100,
+            # 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)": 100.0e00,
+            "patch_icpp(1)%vel(2)": 100.0e00,
+            "patch_icpp(1)%pres": 1.0e05,
+            "patch_icpp(1)%alpha_rho(1)": 1000.0e00,
+            "patch_icpp(1)%alpha_rho(2)": 1.0,
+            "patch_icpp(1)%alpha(1)": 1.0e-12,
+            "patch_icpp(1)%alpha(2)": 1.0 - 1.0e-12,
+            # Patch 2: Rotated rectangle to transport
+            "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)": 100.0e00,
+            "patch_icpp(2)%vel(2)": 100.0e00,
+            "patch_icpp(2)%pres": 1.0e05,
+            "patch_icpp(2)%alpha_rho(1)": 1.0,
+            "patch_icpp(2)%alpha_rho(2)": 1.0,
+            "patch_icpp(2)%alpha(1)": 0,
+            "patch_icpp(2)%alpha(2)": 1.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (2.35e00 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 2.35e00 * 1.0e09 / (2.35e00 - 1.0e00),
+            "fluid_pp(2)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
+            "fluid_pp(2)%pi_inf": 0.0e00,
+        }
+    )
+)

src/common/m_derived_types.fpp

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

src/common/m_patch_geometries.fpp

@@ -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
 
@@ -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
+        real(wp)                                  :: angle
+
+        ! 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

src/pre_process/m_global_parameters.fpp

@@ -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
             patch_icpp(i)%radii = dflt_real
             patch_icpp(i)%alter_patch = .false.
             patch_icpp(i)%alter_patch(0) = .true.

src/pre_process/m_icpp_patches.fpp

@@ -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
@@ -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
@@ -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)
@@ -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()
 
@@ -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
+                ! 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)
 

src/simulation/m_ib_patches.fpp

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

toolchain/mfc/params/definitions.py

@@ -820,6 +820,7 @@ def _load():
         for j in range(1, 4):
             _r(f"{px}radii({j})", REAL)
             _r(f"{px}normal({j})", REAL)
+            _r(f"{px}angles({j})", REAL)
             _r(f"{px}vel({j})", A_REAL, math=r"\f$u_" + str(j) + r"\f$")
         for f in range(1, NF + 1):
             _r(f"{px}alpha({f})", A_REAL, math=r"\f$\alpha_" + str(f) + r"\f$")