Fix reflective BC for stresses

src/common/m_boundary_conditions.fpp

@@ -354,6 +354,13 @@
                                     q_prim_vf(i)%sf(j - 1, k, l)
                             end do
 
+                            if (elasticity) then
+                                do i = 1, shear_BC_flip_num
+                                    q_prim_vf(shear_BC_flip_indices(1, i))%sf(-j, k, l) = &
+                                        -q_prim_vf(shear_BC_flip_indices(1, i))%sf(j - 1, k, l)
+                                end do
+                            end if
+
                             if (hyperelasticity) then
                                 q_prim_vf(xibeg)%sf(-j, k, l) = &
                                     -q_prim_vf(xibeg)%sf(j - 1, k, l)
@@ -404,6 +411,13 @@
                                     q_prim_vf(i)%sf(m - (j - 1), k, l)
                             end do
 
+                            if (elasticity) then
+                                do i = 1, shear_BC_flip_num
+                                    q_prim_vf(shear_BC_flip_indices(1, i))%sf(m + j, k, l) = &
+                                        -q_prim_vf(shear_BC_flip_indices(1, i))%sf(m - (j - 1), k, l)
+                                end do
+                            end if
+
                             if (hyperelasticity) then
                                 q_prim_vf(xibeg)%sf(m + j, k, l) = &
                                     -q_prim_vf(xibeg)%sf(m - (j - 1), k, l)
@@ -457,6 +471,13 @@
                                     q_prim_vf(i)%sf(l, j - 1, k)
                             end do
 
+                            if (elasticity) then
+                                do i = 1, shear_BC_flip_num
+                                    q_prim_vf(shear_BC_flip_indices(2, i))%sf(l, -j, k) = &
+                                        -q_prim_vf(shear_BC_flip_indices(2, i))%sf(l, j - 1, k)
+                                end do
+                            end if
+
                             if (hyperelasticity) then
                                 q_prim_vf(xibeg + 1)%sf(l, -j, k) = &
                                     -q_prim_vf(xibeg + 1)%sf(l, j - 1, k)
@@ -504,6 +525,13 @@
                                     q_prim_vf(i)%sf(l, n - (j - 1), k)
                             end do
 
+                            if (elasticity) then
+                                do i = 1, shear_BC_flip_num
+                                    q_prim_vf(shear_BC_flip_indices(2, i))%sf(l, n + j, k) = &
+                                        -q_prim_vf(shear_BC_flip_indices(2, i))%sf(l, n - (j - 1), k)
+                                end do
+                            end if
+
                             if (hyperelasticity) then
                                 q_prim_vf(xibeg + 1)%sf(l, n + j, k) = &
                                     -q_prim_vf(xibeg + 1)%sf(l, n - (j - 1), k)
@@ -556,6 +584,13 @@
                                     q_prim_vf(i)%sf(k, l, j - 1)
                             end do
 
+                            if (elasticity) then
+                                do i = 1, shear_BC_flip_num
+                                    q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, l, -j) = &
+                                        -q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, l, j - 1)
+                                end do
+                            end if
+
                             if (hyperelasticity) then
                                 q_prim_vf(xiend)%sf(k, l, -j) = &
                                     -q_prim_vf(xiend)%sf(k, l, j - 1)
@@ -603,6 +638,13 @@
                                     q_prim_vf(i)%sf(k, l, p - (j - 1))
                             end do
 
+                            if (elasticity) then
+                                do i = 1, shear_BC_flip_num
+                                    q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, l, p + j) = &
+                                        -q_prim_vf(shear_BC_flip_indices(3, i))%sf(k, l, p - (j - 1))
+                                end do
+                            end if
+
                             if (hyperelasticity) then
                                 q_prim_vf(xiend)%sf(k, l, p + j) = &
                                     -q_prim_vf(xiend)%sf(k, l, p - (j - 1))

src/post_process/m_global_parameters.fpp

@@ -150,6 +150,15 @@
     type(int_bounds_info) :: bc_x, bc_y, bc_z
     !> @}
 
+    integer :: shear_num !! Number of shear stress components
+    integer, dimension(3) :: shear_indices !<
+    !! Indices of the stress components that represent shear stress
+    integer :: shear_BC_flip_num !<
+    !! Number of shear stress components to reflect for boundary conditions
+    integer, dimension(3, 2) :: shear_BC_flip_indices !<
+    !! Indices of shear stress components to reflect for boundary conditions.
+    !! Size: (1:3, 1:shear_BC_flip_num) for (x/y/z, [indices])
+
     logical :: parallel_io    !< Format of the data files
     logical :: sim_data
     logical :: file_per_process !< output format
@@ -585,29 +594,6 @@
                 sys_size = B_idx%end
             end if
 
-            if (hypoelasticity .or. hyperelasticity) then
-                elasticity = .true.
-                stress_idx%beg = sys_size + 1
-                stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
-                ! number of distinct stresses is 1 in 1D, 3 in 2D, 6 in 3D
-                sys_size = stress_idx%end
-            end if
-
-            if (hyperelasticity) then
-                xi_idx%beg = sys_size + 1
-                xi_idx%end = sys_size + num_dims
-                ! adding three more equations for the \xi field and the elastic energy
-                sys_size = xi_idx%end + 1
-                ! number of entries in the symmetric btensor plus the jacobian
-                b_size = (num_dims*(num_dims + 1))/2 + 1
-                tensor_size = num_dims**2 + 1
-            end if
-
-            if (surface_tension) then
-                c_idx = sys_size + 1
-                sys_size = c_idx
-            end if
-
             ! Volume Fraction Model (6-equation model)
         else if (model_eqns == 3) then
 
@@ -626,29 +612,6 @@
             sys_size = internalEnergies_idx%end
             alf_idx = 1 ! dummy, cannot actually have a void fraction
 
-            if (hypoelasticity .or. hyperelasticity) then
-                elasticity = .true.
-                stress_idx%beg = sys_size + 1
-                stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
-                ! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
-                sys_size = stress_idx%end
-            end if
-
-            if (hyperelasticity) then
-                xi_idx%beg = sys_size + 1
-                xi_idx%end = sys_size + num_dims
-                ! adding three more equations for the \xi field and the elastic energy
-                sys_size = xi_idx%end + 1
-                ! number of entries in the symmetric btensor plus the jacobian
-                b_size = (num_dims*(num_dims + 1))/2 + 1
-                tensor_size = num_dims**2 + 1
-            end if
-
-            if (surface_tension) then
-                c_idx = sys_size + 1
-                sys_size = c_idx
-            end if
-
         else if (model_eqns == 4) then
             cont_idx%beg = 1 ! one continuity equation
             cont_idx%end = 1 !num_fluids
@@ -705,6 +668,54 @@
             end if
         end if
 
+        if (model_eqns == 2 .or. model_eqns == 3) then
+
+            if (hypoelasticity .or. hyperelasticity) then
+                elasticity = .true.
+                stress_idx%beg = sys_size + 1
+                stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
+                ! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
+                sys_size = stress_idx%end
+
+                ! shear stress index is 2 for 2D and 2,4,5 for 3D
+                if (num_dims == 1) then
+                    shear_num = 0
+                else if (num_dims == 2) then
+                    shear_num = 1
+                    shear_indices(1) = stress_idx%beg - 1 + 2
+                    shear_BC_flip_num = 1
+                    shear_BC_flip_indices(1:2, 1) = shear_indices(1)
+                    ! Both x-dir and y-dir: flip tau_xy only
+                else if (num_dims == 3) then
+                    shear_num = 3
+                    shear_indices(1:3) = stress_idx%beg - 1 + (/2, 4, 5/)
+                    shear_BC_flip_num = 2
+                    shear_BC_flip_indices(1, 1:2) = shear_indices((/1, 2/))
+                    shear_BC_flip_indices(2, 1:2) = shear_indices((/1, 3/))
+                    shear_BC_flip_indices(3, 1:2) = shear_indices((/2, 3/))
+                    ! x-dir: flip tau_xy and tau_xz
+                    ! y-dir: flip tau_xy and tau_yz
+                    ! z-dir: flip tau_xz and tau_yz
+                end if
+            end if
+
+            if (hyperelasticity) then
+                xi_idx%beg = sys_size + 1
+                xi_idx%end = sys_size + num_dims
+                ! adding three more equations for the \xi field and the elastic energy
+                sys_size = xi_idx%end + 1
+                ! number of entries in the symmetric btensor plus the jacobian
+                b_size = (num_dims*(num_dims + 1))/2 + 1
+                tensor_size = num_dims**2 + 1
+            end if
+
+            if (surface_tension) then
+                c_idx = sys_size + 1
+                sys_size = c_idx
+            end if
+
+        end if
+
         if (chemistry) then
             species_idx%beg = sys_size + 1
             species_idx%end = sys_size + num_species

src/pre_process/m_global_parameters.fpp

@@ -124,6 +124,15 @@
     type(int_bounds_info) :: bc_x, bc_y, bc_z !<
     !! Boundary conditions in the x-, y- and z-coordinate directions
 
+    integer :: shear_num !! Number of shear stress components
+    integer, dimension(3) :: shear_indices !<
+    !! Indices of the stress components that represent shear stress
+    integer :: shear_BC_flip_num !<
+    !! Number of shear stress components to reflect for boundary conditions
+    integer, dimension(3, 2) :: shear_BC_flip_indices !<
+    !! Indices of shear stress components to reflect for boundary conditions.
+    !! Size: (1:3, 1:shear_BC_flip_num) for (x/y/z, [indices])
+
     logical :: parallel_io !< Format of the data files
     logical :: file_per_process !< type of data output
     integer :: precision !< Precision of output files
@@ -677,29 +686,6 @@
                 sys_size = B_idx%end
             end if
 
-            if (hypoelasticity .or. hyperelasticity) then
-                elasticity = .true.
-                stress_idx%beg = sys_size + 1
-                stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
-                ! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
-                sys_size = stress_idx%end
-            end if
-
-            if (hyperelasticity) then
-                ! number of entries in the symmetric btensor plus the jacobian
-                b_size = (num_dims*(num_dims + 1))/2 + 1
-                tensor_size = num_dims**2 + 1
-                xi_idx%beg = sys_size + 1
-                xi_idx%end = sys_size + num_dims
-                ! adding three more equations for the \xi field and the elastic energy
-                sys_size = xi_idx%end + 1
-            end if
-
-            if (surface_tension) then
-                c_idx = sys_size + 1
-                sys_size = c_idx
-            end if
-
             ! Volume Fraction Model (6-equation model)
         else if (model_eqns == 3) then
 
@@ -717,29 +703,6 @@
             internalEnergies_idx%end = adv_idx%end + num_fluids
             sys_size = internalEnergies_idx%end
 
-            if (hypoelasticity .or. hyperelasticity) then
-                elasticity = .true.
-                stress_idx%beg = sys_size + 1
-                stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
-                ! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
-                sys_size = stress_idx%end
-            end if
-
-            if (hyperelasticity) then
-                ! number of entries in the symmetric btensor plus the jacobian
-                b_size = (num_dims*(num_dims + 1))/2 + 1
-                tensor_size = num_dims**2 + 1
-                xi_idx%beg = sys_size + 1
-                xi_idx%end = sys_size + num_dims
-                ! adding three more equations for the \xi field and the elastic energy
-                sys_size = xi_idx%end + 1
-            end if
-
-            if (surface_tension) then
-                c_idx = sys_size + 1
-                sys_size = c_idx
-            end if
-
         else if (model_eqns == 4) then
             ! 4 equation model with subgrid bubbles_euler
             cont_idx%beg = 1 ! one continuity equation
@@ -798,6 +761,54 @@
             end if
         end if
 
+        if (model_eqns == 2 .or. model_eqns == 3) then
+
+            if (hypoelasticity .or. hyperelasticity) then
+                elasticity = .true.
+                stress_idx%beg = sys_size + 1
+                stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
+                ! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
+                sys_size = stress_idx%end
+
+                ! shear stress index is 2 for 2D and 2,4,5 for 3D
+                if (num_dims == 1) then
+                    shear_num = 0
+                else if (num_dims == 2) then
+                    shear_num = 1
+                    shear_indices(1) = stress_idx%beg - 1 + 2
+                    shear_BC_flip_num = 1
+                    shear_BC_flip_indices(1:2, 1) = shear_indices(1)
+                    ! Both x-dir and y-dir: flip tau_xy only
+                else if (num_dims == 3) then
+                    shear_num = 3
+                    shear_indices(1:3) = stress_idx%beg - 1 + (/2, 4, 5/)
+                    shear_BC_flip_num = 2
+                    shear_BC_flip_indices(1, 1:2) = shear_indices((/1, 2/))
+                    shear_BC_flip_indices(2, 1:2) = shear_indices((/1, 3/))
+                    shear_BC_flip_indices(3, 1:2) = shear_indices((/2, 3/))
+                    ! x-dir: flip tau_xy and tau_xz
+                    ! y-dir: flip tau_xy and tau_yz
+                    ! z-dir: flip tau_xz and tau_yz
+                end if
+            end if
+
+            if (hyperelasticity) then
+                ! number of entries in the symmetric btensor plus the jacobian
+                b_size = (num_dims*(num_dims + 1))/2 + 1
+                tensor_size = num_dims**2 + 1
+                xi_idx%beg = sys_size + 1
+                xi_idx%end = sys_size + num_dims
+                ! adding three more equations for the \xi field and the elastic energy
+                sys_size = xi_idx%end + 1
+            end if
+
+            if (surface_tension) then
+                c_idx = sys_size + 1
+                sys_size = c_idx
+            end if
+
+        end if
+
         if (chemistry) then
             species_idx%beg = sys_size + 1
             species_idx%end = sys_size + num_species