feature: integrate rotations into ICPP patches (#1594)

examples/2D_rotated_fluid_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": 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,
+        }
+    )
+)

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

src/common/m_model.fpp

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

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,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
+
+        ! 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

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