Fix STL+IBM test suite and refactor

Author: haochey

examples/2D_IBM_STL/case.py

@@ -6,10 +6,10 @@
 # Dynamic Viscosity
 Mu1 = 0.0000184
 #Mu2 = 0.01
-rho1 = 1.19 #rho1 = 0.2199
+rho1 = 1.19 #0.2199
 gam_a = 1.4
 # Patch Design
-D = 0.01
+D = 0.1
 
 # Configuring case dictionary
 print(json.dumps({
@@ -22,10 +22,10 @@
     'x_domain%end'                 : 6*D,
     'y_domain%beg'                 : -3*D,
     'y_domain%end'                 : 3*D,
-    'm'                            : 199,
-    'n'                            : 99,
+    'm'                            : 60,
+    'n'                            : 30,
     'p'                            : 0,
-    'dt'                           : 1.0E-12,
+    'dt'                           : 1.0E-9,
     't_step_start'                 : 0,
     't_step_stop'                  : 3000,
     't_step_save'                  : 30,
@@ -47,7 +47,7 @@
     'mp_weno'                      : 'F',
     'riemann_solver'               : 2,
     'wave_speeds'                  : 1,
-    'viscous'                      : 'T',
+    'viscous'                      : 'F',
     'bc_x%beg'                     : -3,
     'bc_x%end'                     : -3,
     'bc_y%beg'                     : -3,
@@ -71,26 +71,24 @@
     'patch_icpp(1)%y_centroid'     : 0,
     'patch_icpp(1)%length_x'       : 1000*D,
     'patch_icpp(1)%length_y'       : 1000*D,
-    'patch_icpp(1)%vel(1)'         : 10, #527.2E+00,
+    'patch_icpp(1)%vel(1)'         : 1,
     'patch_icpp(1)%vel(2)'         : 0.0E+00,
-    'patch_icpp(1)%pres'           : 100000, #10918.2549,
+    'patch_icpp(1)%pres'           : 100000,
     'patch_icpp(1)%alpha_rho(1)'   : (1.0)*rho1,
     'patch_icpp(1)%alpha(1)'       : 1.0,
     # ==========================================================================
     'patch_ib(1)%geometry'              : 5,
     'patch_ib(1)%model%filepath'        : 'Circle_IBM.stl',
-    'patch_ib(1)%model%translate(1)'    : -0.05 + 0.0424444350,
-    'patch_ib(1)%model%translate(2)'    : -0.05 + 0.0445000000,
-    'patch_ib(1)%model%spc'             : 2000,
+    'patch_ib(1)%model%translate(1)'    : -0.05,
+    'patch_ib(1)%model%translate(2)'    : -0.05,
+    'patch_ib(1)%model%spc'             : 200,
     'patch_ib(1)%model%threshold'       : 0.99,
     'patch_ib(1)%slip'                  : 'F',
     # # ========================================================================
 
     # Fluids Physical Parameters ===============================================
     'fluid_pp(1)%gamma'            : 1.E+00/(gam_a-1.E+00),
     'fluid_pp(1)%pi_inf'           : 0,
-    'fluid_pp(1)%Re(1)'            : 5000,
-    #'fluid_pp(1)%Re(1)'            : 7535533.2,
     # ==========================================================================
 }))
 

examples/3D_IBM_STL/case.py

@@ -3,8 +3,8 @@
 
 Mu = 1.84E-05
 gam_a = 1.4
-rho1 = 0.2199
-D = 0.02
+rho1 = 1.19
+D = 2
 
 # Configuring case dictionary
 print(json.dumps({
@@ -13,8 +13,7 @@
     # ==========================================================================
 
     # Computational Domain Parameters ==========================================
-    # x direction
-    'x_domain%beg'                 : -4.0*D,
+    'x_domain%beg'                 : -4*D,
     'x_domain%end'                 : 4.0*D,
     # y direction
     'y_domain%beg'                 : -2*D,
@@ -23,10 +22,10 @@
     'z_domain%beg'                 : -2*D,
     'z_domain%end'                 : 2*D,
     'cyl_coord'                    : 'F',
-    'm'                            : 99,
-    'n'                            : 49,
-    'p'                            : 49,
-    'dt'                           : 1.0E-12,
+    'm'                            : 59,
+    'n'                            : 29,
+    'p'                            : 29,
+    'dt'                           : 1.0E-9,
     't_step_start'                 : 0,
     't_step_stop'                  : 1000,
     't_step_save'                  : 10,
@@ -42,15 +41,15 @@
     'time_stepper'                 : 3,
     'weno_order'                   : 5,
     'weno_eps'                     : 1.E-16,
-    'weno_Re_flux'                 : 'T',
+    'weno_Re_flux'                 : 'F',
     'weno_avg'                     : 'T',
     'avg_state'                    : 2,
     'mapped_weno'                  : 'T',
     'null_weights'                 : 'F',
     'mp_weno'                      : 'F',
     'riemann_solver'               : 2,
     'wave_speeds'                  : 1,
-    'viscous'                      :'T',
+    'viscous'                      :'F',
     'bc_x%beg'                     : -3,
     'bc_x%end'                     : -3,
     'bc_y%beg'                     : -3,
@@ -81,28 +80,27 @@
     'patch_icpp(1)%length_x'       : 100*D,
     'patch_icpp(1)%length_y'       : 50*D,
     'patch_icpp(1)%length_z'       : 50*D,
-    'patch_icpp(1)%vel(1)'         : 527.2E+00,
+    'patch_icpp(1)%vel(1)'         : 1,
     'patch_icpp(1)%vel(2)'         : 0.0E+00,
     'patch_icpp(1)%vel(3)'         : 0.0E+00,
-    'patch_icpp(1)%pres'           : 10918.2549,
+    'patch_icpp(1)%pres'           : 100000,
     'patch_icpp(1)%alpha_rho(1)'   : (1.0)*rho1,
     'patch_icpp(1)%alpha(1)'       : 1.E+00,
     # # ========================================================================
 
     # Patch: Model Immersed Boundary ===========================================
     'patch_ib(1)%geometry'                     : 12,
     'patch_ib(1)%model%filepath'               : 'Sphere_IBM.stl',
-    'patch_ib(1)%model%translate(1)'           : 0 - 0.02056,
-    'patch_ib(1)%model%translate(2)'           : 0 - 0.01,
-    'patch_ib(1)%model%translate(3)'           : 0 - 0.01,
-    'patch_ib(1)%model%spc'                    : 200,
-    'patch_ib(1)%model%threshold'              : 0.99,
+    'patch_ib(1)%model%translate(1)'           : 0,
+    'patch_ib(1)%model%translate(2)'           : -0.05,
+    'patch_ib(1)%model%translate(3)'           : -0.05,
+    'patch_ib(1)%model%spc'                    : 100,
+    'patch_ib(1)%model%threshold'              : 0.95,
     'patch_ib(1)%slip'                         : 'F',
     # ==========================================================================
 
     # Fluids Physical Parameters ===============================================
     'fluid_pp(1)%gamma'            : 1.E+00/(gam_a-1.E+00),
     'fluid_pp(1)%pi_inf'           : 0,
-    'fluid_pp(1)%Re(1)'            : 7535533.2,
     # ==========================================================================
 }))

src/common/m_constants.fpp

@@ -32,6 +32,8 @@ module m_constants
     integer, parameter :: Ifactor_2D = 100 !< Multiple factor of the ratio (edge to cell width) for interpolation along edges for 2D models
     integer, parameter :: Ifactor_3D = 10 !< Multiple factor of the ratio (edge to cell width) for interpolation along edges for 3D models
     integer, parameter :: Ifactor_bary_3D = 30 !< Multiple factor of the ratio (triangle area to cell face area) for interpolation on triangle facets for 3D models
+    integer, parameter :: num_ray = 20 !< Default number of rays traced per cell
+    real(kind(0d0)), parameter :: ray_tracing_threshold = 0.9d0 !< Threshold above which the cell is marked as the model patch
     real(kind(0d0)), parameter :: threshold_vector_zero = 1d-10 !< Threshold to treat the component of a vector to be zero
     real(kind(0d0)), parameter :: threshold_edge_zero = 1d-8 !< Threshold to treat two edges to be overlapped
     real(kind(0d0)), parameter :: threshold_bary = 1d-1 !< Threshold to interpolate a barycentric facet

src/pre_process/m_check_ib_patches.fpp

@@ -269,7 +269,7 @@ contains
 
         call s_int_to_str(patch_id, iStr)
 
-        @:PROHIBIT(patch_ib(patch_id)%model%filepath == ' ', &
+        @:PROHIBIT(patch_ib(patch_id)%model%filepath == dflt_char, &
             'Empty model file path for patch '//trim(iStr))
 
         @:PROHIBIT(patch_ib(patch_id)%model%scale(1) <= 0d0 &

src/pre_process/m_global_parameters.fpp

@@ -342,9 +342,9 @@ contains
             patch_icpp(i)%geometry = dflt_int
             patch_icpp(i)%model%scale(:) = 1d0
             patch_icpp(i)%model%translate(:) = 0d0
-            patch_icpp(i)%model%filepath(:) = ' '
-            patch_icpp(i)%model%spc = 10
-            patch_icpp(i)%model%threshold = 0.9d0
+            patch_icpp(i)%model%filepath(:) = dflt_char
+            patch_icpp(i)%model%spc = num_ray
+            patch_icpp(i)%model%threshold = ray_tracing_threshold
             patch_icpp(i)%x_centroid = dflt_real
             patch_icpp(i)%y_centroid = dflt_real
             patch_icpp(i)%z_centroid = dflt_real
@@ -450,9 +450,9 @@ contains
             patch_ib(i)%model%scale(:) = 1d0
             patch_ib(i)%model%translate(:) = 0d0
             patch_ib(i)%model%rotate(:) = 0d0
-            patch_ib(i)%model%filepath(:) = ' '
-            patch_ib(i)%model%spc = 20
-            patch_ib(i)%model%threshold = 0.9d0
+            patch_ib(i)%model%filepath(:) = dflt_char
+            patch_ib(i)%model%spc = num_ray
+            patch_ib(i)%model%threshold = ray_tracing_threshold
         end do
 
         ! Fluids physical parameters

src/pre_process/m_initial_condition.fpp

@@ -229,11 +229,11 @@ contains
                 !> @{
                 ! Circular patch
                 if (patch_icpp(i)%geometry == 2) then
-                    call s_circle(i, patch_id_fp, q_prim_vf, .false.)
+                    call s_circle(i, patch_id_fp, q_prim_vf)
 
                     ! Rectangular patch
                 elseif (patch_icpp(i)%geometry == 3) then
-                    call s_rectangle(i, patch_id_fp, q_prim_vf, .false.)
+                    call s_rectangle(i, patch_id_fp, q_prim_vf)
 
                     ! Swept line patch
                 elseif (patch_icpp(i)%geometry == 4) then
@@ -266,7 +266,7 @@ contains
 
                     ! STL patch
                 elseif (patch_icpp(i)%geometry == 21) then
-                    call s_model(i, patch_id_fp, q_prim_vf, .false.)
+                    call s_model(i, patch_id_fp, q_prim_vf)
 
                 end if
                 !> @}

src/pre_process/m_model.fpp

@@ -582,18 +582,18 @@ contains
         ! Collect all edges of all triangles and store them
         do i = 1, model%ntrs
             ! First edge (v1, v2)
-            edge(1, :) = model%trs(i)%v(1, 1:2)
-            edge(2, :) = model%trs(i)%v(2, 1:2)
+            edge(1, 1:2) = model%trs(i)%v(1, 1:2)
+            edge(2, 1:2) = model%trs(i)%v(2, 1:2)
             call f_register_edge(temp_boundary_v, edge, edge_index, edge_count)
 
             ! Second edge (v2, v3)
-            edge(1, :) = model%trs(i)%v(2, 1:2)
-            edge(2, :) = model%trs(i)%v(3, 1:2)
+            edge(1, 1:2) = model%trs(i)%v(2, 1:2)
+            edge(2, 1:2) = model%trs(i)%v(3, 1:2)
             call f_register_edge(temp_boundary_v, edge, edge_index, edge_count)
 
             ! Third edge (v3, v1)
-            edge(1, :) = model%trs(i)%v(3, 1:2)
-            edge(2, :) = model%trs(i)%v(1, 1:2)
+            edge(1, 1:2) = model%trs(i)%v(3, 1:2)
+            edge(2, 1:2) = model%trs(i)%v(1, 1:2)
             call f_register_edge(temp_boundary_v, edge, edge_index, edge_count)
         end do
 
@@ -635,8 +635,8 @@ contains
         do i = 1, edge_count
             if (edge_occurrence(i) == 0) then
                 store_index = store_index + 1
-                boundary_v(store_index, 1, :) = temp_boundary_v(i, 1, :)
-                boundary_v(store_index, 2, :) = temp_boundary_v(i, 2, :)
+                boundary_v(store_index, 1, 1:2) = temp_boundary_v(i, 1, 1:2)
+                boundary_v(store_index, 2, 1:2) = temp_boundary_v(i, 2, 1:2)
             end if
         end do
 
@@ -680,8 +680,8 @@ contains
 
         ! Increment edge index and store the edge
         edge_index = edge_index + 1
-        temp_boundary_v(edge_index, 1, :) = edge(1, :)
-        temp_boundary_v(edge_index, 2, :) = edge(2, :)
+        temp_boundary_v(edge_index, 1, 1:2) = edge(1, 1:2)
+        temp_boundary_v(edge_index, 2, 1:2) = edge(2, 1:2)
 
     end subroutine f_register_edge
 
@@ -693,24 +693,19 @@ contains
     subroutine f_check_interpolation_2D(boundary_v, boundary_edge_count, spacing, interpolate)
         logical, intent(out) :: interpolate !< Logical indicator of interpolation
         integer, intent(in) :: boundary_edge_count !< Number of boundary edges
-        real(kind(0d0)), optional, intent(in), dimension(1:boundary_edge_count, 1:3, 1:2) :: boundary_v
+        real(kind(0d0)), intent(in), dimension(1:boundary_edge_count, 1:3, 1:2) :: boundary_v
         t_vec3, intent(in) :: spacing
 
         real(kind(0d0)) :: l1, cell_width !< Length of each boundary edge and cell width
-        real(kind(0d0)), dimension(1:2, 1:2) :: edge_v !< Boundary points of each boundary edge
         integer :: j !< Boundary edge index iterator
 
         cell_width = minval(spacing(1:2))
         interpolate = .false.
 
         do j = 1, boundary_edge_count
-            edge_v(1, 1) = boundary_v(j, 1, 1)
-            edge_v(2, 1) = boundary_v(j, 2, 1)
-            edge_v(1, 2) = boundary_v(j, 1, 2)
-            edge_v(2, 2) = boundary_v(j, 2, 2)
 
-            l1 = dsqrt((edge_v(2, 1) - edge_v(1, 1))**2 + &
-                       (edge_v(2, 2) - edge_v(1, 2))**2)
+            l1 = dsqrt((boundary_v(j, 2, 1) - boundary_v(j, 1, 1))**2 + &
+                       (boundary_v(j, 2, 2) - boundary_v(j, 1, 2))**2)
 
             if ((l1 > cell_width)) then
                 interpolate = .true.

src/pre_process/m_patches.fpp

@@ -2012,7 +2012,7 @@ contains
         if (p > 0) then
             call f_check_interpolation_3D(model, (/dx, dy, dz/), interpolate)
         else
-            call f_check_interpolation_2D(boundary_v, boundary_vertex_count, (/dx, dy, dz/), interpolate)
+            call f_check_interpolation_2D(boundary_v, boundary_edge_count, (/dx, dy, dz/), interpolate)
         end if
 
         ! Show the number of edges and boundary edges in 2D STL models