Viscous stress and ellipse ib

examples/2D_ibm_ellipse/case.py

@@ -0,0 +1,97 @@
+import json
+import math
+
+Mu = 1.84e-05
+gam_a = 1.4
+
+# Configuring case dictionary
+print(
+    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": 250,
+            "n": 250,
+            "p": 0,
+            "dt": 0.5e-5,
+            "t_step_start": 0,
+            "t_step_stop": 1000,  # 3000
+            "t_step_save": 10,  # 10
+            # 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": 6,
+            "patch_ib(1)%x_centroid": 1.5e-03,
+            "patch_ib(1)%y_centroid": 3.0e-03,
+            "patch_ib(1)%length_x": 0.4e-03,
+            "patch_ib(1)%length_y": 0.2e-03,
+            "patch_ib(1)%slip": "F",
+            # 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,
+        }
+    )
+)

examples/2D_mibm_shock_cylinder/case.py

@@ -0,0 +1,138 @@
+import json
+import math
+
+# This case is a recreation of the case from "Moving overlapping grids with adaptive mesh refinement for high-speed reactive and non-reactive flow"
+# by William D. Henshaw and Donald W. Schwendeman
+
+# fluid parameters
+gam_a = 1.4
+
+# domain size and speed
+mach_number = 1.5
+pre_shock_pressure = 1
+pre_shock_density = 1.4
+pre_shock_speed = 0.0
+post_shock_pressure = 2.4583
+post_shock_density = 2.6069
+post_shock_speed = 0.6944
+
+domain_size = 4.0
+wave_front = -1.5
+
+total_time = 1.5
+num_time_steps = 2000
+dt = float(total_time / num_time_steps)
+num_saves = 100
+steps_to_save = int(num_time_steps / num_saves)
+
+# Configuring case dictionary
+print(
+    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": -domain_size * 0.5,
+            "x_domain%end": domain_size * 0.5,
+            # r direction
+            "y_domain%beg": -domain_size * 0.5,
+            "y_domain%end": domain_size * 0.5,
+            "cyl_coord": "F",
+            "m": 1000,
+            "n": 1000,
+            "p": 0,
+            "dt": dt,
+            "t_step_start": 0,
+            "t_step_stop": num_time_steps,  # 10000,
+            "t_step_save": steps_to_save,
+            # Simulation Algorithm Parameters
+            # Only one patches are necessary, the air tube
+            "num_patches": 2,
+            # 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": -17,
+            "bc_x%end": -8,
+            "bc_y%beg": -15,
+            "bc_y%end": -15,
+            # 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(2)%geometry": 3,
+            # patch locations
+            "patch_icpp(1)%x_centroid": 0.5 * wave_front + 0.25 * domain_size,
+            "patch_icpp(1)%y_centroid": 0.0,
+            "patch_icpp(1)%length_x": 0.5 * domain_size - wave_front,
+            "patch_icpp(1)%length_y": domain_size,
+            "patch_icpp(2)%x_centroid": 0.5 * wave_front - 0.25 * domain_size,
+            "patch_icpp(2)%y_centroid": 0.0,
+            "patch_icpp(2)%length_x": 0.5 * domain_size + wave_front,
+            "patch_icpp(2)%length_y": domain_size,
+            # Specify the patch primitive variables
+            "patch_icpp(1)%vel(1)": pre_shock_speed,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": pre_shock_pressure,
+            "patch_icpp(1)%alpha_rho(1)": pre_shock_density,
+            "patch_icpp(1)%alpha(1)": 1.0e00,
+            "patch_icpp(2)%vel(1)": post_shock_speed,
+            "patch_icpp(2)%vel(2)": 0.0,
+            "patch_icpp(2)%pres": post_shock_pressure,
+            "patch_icpp(2)%alpha_rho(1)": post_shock_density,
+            "patch_icpp(2)%alpha(1)": 1.0e00,
+            # Patch: Cylinder Immersed Boundary
+            "patch_ib(1)%geometry": 2,
+            "patch_ib(1)%x_centroid": -0.5,
+            "patch_ib(1)%y_centroid": 0.0,
+            "patch_ib(1)%radius": 0.5,
+            "patch_ib(1)%slip": "T",
+            "patch_ib(1)%moving_ibm": 2,
+            "patch_ib(1)%vel(1)": 0,
+            "patch_ib(1)%vel(2)": 0,
+            "patch_ib(1)%vel(3)": 0,
+            "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": 0.5,  # z-axis rotation
+            # 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,
+        }
+    )
+)

src/common/m_compute_levelset.fpp

@@ -23,7 +23,8 @@ module m_compute_levelset
  s_3D_airfoil_levelset, &
  s_rectangle_levelset, &
  s_cuboid_levelset, &
- s_sphere_levelset
+ s_sphere_levelset, &
+ s_ellipse_levelset
 
 contains
 
@@ -336,6 +337,62 @@ contains
 
     end subroutine s_rectangle_levelset
 
+    subroutine s_ellipse_levelset(ib_patch_id, levelset, levelset_norm)
+
+        type(levelset_field), intent(INOUT), optional :: levelset
+        type(levelset_norm_field), intent(INOUT), optional :: levelset_norm
+
+        integer, intent(in) :: ib_patch_id
+        real(wp) :: ellipse_coeffs(2) ! a and b in the ellipse equation
+        real(wp) :: quadratic_coeffs(3) ! A, B, C in the quadratic equation to compute levelset
+
+        real(wp) :: length_x, length_y
+        real(wp), dimension(1:3) :: xy_local, normal_vector !< x and y coordinates in local IB frame
+        real(wp), dimension(2) :: center !< x and y coordinates in local IB frame
+        real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
+
+        integer :: i, j, k !< Loop index variables
+        integer :: idx !< Shortest path direction indicator
+
+        length_x = patch_ib(ib_patch_id)%length_x
+        length_y = patch_ib(ib_patch_id)%length_y
+        center(1) = patch_ib(ib_patch_id)%x_centroid
+        center(2) = patch_ib(ib_patch_id)%y_centroid
+        inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
+        rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)
+
+        ellipse_coeffs(1) = 0.5_wp*length_x
+        ellipse_coeffs(2) = 0.5_wp*length_y
+
+        $:GPU_PARALLEL_LOOP(private='[i,j,k,idx,quadratic_coeffs,xy_local,normal_vector]', &
+                  & copyin='[ib_patch_id,center,ellipse_coeffs,inverse_rotation,rotation]', collapse=2)
+        do i = 0, m
+            do j = 0, n
+                xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
+                xy_local = matmul(inverse_rotation, xy_local)
+
+                ! we will get NaNs in the levelset if we compute this outside the ellipse
+                if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
+
+                    normal_vector = xy_local
+                    normal_vector(2) = normal_vector(2)*(ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp ! get the normal direction via the coordinate transformation method
+                    normal_vector = normal_vector/sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector
+                    levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, normal_vector) ! save after rotating the vector to the global frame
+
+                    ! use the normal vector to set up the quadratic equation for the levelset, using A, B, and C in indices 1, 2, and 3
+                    quadratic_coeffs(1) = (normal_vector(1)/ellipse_coeffs(1))**2 + (normal_vector(2)/ellipse_coeffs(2))**2
+                    quadratic_coeffs(2) = 2._wp*((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2)/(ellipse_coeffs(2)**2)))
+                    quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp
+
+                    ! compute the levelset with the quadratic equation [ -B + sqrt(B^2 - 4AC) ] / 2A
+                    levelset%sf(i, j, 0, ib_patch_id) = -0.5_wp*(-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp*quadratic_coeffs(1)*quadratic_coeffs(3)))/quadratic_coeffs(1)
+                end if
+            end do
+        end do
+        $:END_GPU_PARALLEL_LOOP()
+
+    end subroutine s_ellipse_levelset
+
     subroutine s_cuboid_levelset(ib_patch_id, levelset, levelset_norm)
 
         type(levelset_field), intent(INOUT), optional :: levelset

src/common/m_ib_patches.fpp

@@ -119,6 +119,9 @@ contains
                     ! STL+IBM patch
                 elseif (patch_ib(i)%geometry == 5) then
                     call s_ib_model(i, ib_markers_sf, levelset, levelset_norm)
+                elseif (patch_ib(i)%geometry == 6) then
+                    call s_ib_ellipse(i, ib_markers_sf)
+                    call s_ellipse_levelset(i, levelset, levelset_norm)
                 end if
             end do
             !> @}
@@ -762,6 +765,45 @@ contains
 
     end subroutine s_ib_cylinder
 
+    subroutine s_ib_ellipse(patch_id, ib_markers_sf)
+
+        integer, intent(in) :: patch_id
+        integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
+
+        integer :: i, j, k !< generic loop iterators
+        real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
+        real(wp), dimension(1:2) :: ellipse_coeffs !< a and b in the ellipse coefficients
+        real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame
+        real(wp), dimension(1:3, 1:3) :: inverse_rotation
+
+        ! Transferring the ellipse's centroid and length information
+        center(1) = patch_ib(patch_id)%x_centroid
+        center(2) = patch_ib(patch_id)%y_centroid
+        ellipse_coeffs(1) = 0.5_wp*patch_ib(patch_id)%length_x
+        ellipse_coeffs(2) = 0.5_wp*patch_ib(patch_id)%length_y
+        inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
+
+        ! Checking whether the ellipse covers a particular cell in the
+        ! domain
+        $:GPU_PARALLEL_LOOP(private='[i,j, xy_local]', copy='[ib_markers_sf]',&
+                  & copyin='[patch_id,center,ellipse_coeffs,inverse_rotation,x_cc,y_cc]', collapse=2)
+        do j = 0, n
+            do i = 0, m
+                ! get the x and y coordinates in the local IB frame
+                xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
+                xy_local = matmul(inverse_rotation, xy_local)
+
+                ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1
+                if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
+                    ! Updating the patch identities bookkeeping variable
+                    ib_markers_sf(i, j, 0) = patch_id
+                end if
+            end do
+        end do
+        $:END_GPU_PARALLEL_LOOP()
+
+    end subroutine s_ib_ellipse
+
     !> The STL patch is a 2/3D geometry that is imported from an STL file.
     !! @param patch_id is the patch identifier
     !! @param ib_markers_sf Array to track patch ids