1104 mibm higher order central differencing

examples/2D_backward_facing_step/case.py

@@ -54,6 +54,7 @@
             "bc_y%end": -3,
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             "viscous": "T",
             # Formatted Database Files Structure Parameters
             "format": "silo",

examples/2D_forward_facing_step/case.py

@@ -52,6 +52,7 @@
             "bc_y%end": -2,
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             # Formatted Database Files Structure Parameters
             "format": "silo",
             "precision": "double",

examples/2D_ibm/case.py

@@ -60,6 +60,7 @@
             # Set IB to True and add 1 patch
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             "viscous": "T",
             # Formatted Database Files Structure Parameters
             "format": "silo",

examples/2D_ibm_airfoil/case.py

@@ -64,6 +64,7 @@
             # Set IB to True and add 1 patch
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             # Formatted Database Files Structure Parameters
             "format": "silo",
             "precision": "double",

examples/2D_ibm_cfl_dt/case.py

@@ -62,6 +62,7 @@
             # Set IB to True and add 1 patch
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             "viscous": "T",
             # Formatted Database Files Structure Parameters
             "format": "silo",

examples/2D_ibm_ellipse/case.py

@@ -60,6 +60,7 @@
             # Set IB to True and add 1 patch
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             "viscous": "T",
             # Formatted Database Files Structure Parameters
             "format": "silo",

examples/2D_ibm_multiphase/case.py

@@ -59,6 +59,7 @@
             # Set IB to True and add 1 patch
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             # Formatted Database Files Structure Parameters
             "format": "silo",
             "precision": "double",

examples/2D_ibm_stl_MFCCharacter/case.py

@@ -50,6 +50,7 @@
             "bc_y%end": -3,
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             # Formatted Database Files Structure Parameters
             # Export primitive variables in double precision with parallel
             # I/O to minimize I/O computational time during large simulations

examples/2D_ibm_stl_test/case.py

@@ -51,6 +51,7 @@
             "bc_y%end": -3,
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             # Formatted Database Files Structure Parameters
             # Export primitive variables in double precision with parallel
             # I/O to minimize I/O computational time during large simulations

examples/2D_ibm_stl_wedge/case.py

@@ -51,6 +51,7 @@
             "bc_y%end": -3,
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             # Formatted Database Files Structure Parameters
             # Export primitive variables in double precision with parallel
             # I/O to minimize I/O computational time during large simulations

examples/2D_mibm_cylinder_in_cross_flow/case.py

@@ -52,6 +52,7 @@
             "mp_weno": "T",
             "riemann_solver": "hllc",
             "wave_speeds": "direct",
+            "fd_order": 2,
             # We use ghost-cell
             "bc_x%beg": -3,
             "bc_x%end": -3,

examples/2D_mibm_particle_cloud/case.py

@@ -71,6 +71,7 @@
             "mp_weno": "T",
             "riemann_solver": "hllc",
             "wave_speeds": "direct",
+            "fd_order": 2,
             "bc_x%beg": -17,
             "bc_x%end": -8,
             "bc_y%beg": -15,

examples/2D_mibm_shock_cylinder/case.py

@@ -73,6 +73,7 @@
             "mp_weno": "T",
             "riemann_solver": "hllc",
             "wave_speeds": "direct",
+            "fd_order": 2,
             # We use ghost-cell
             "bc_x%beg": -17,
             "bc_x%end": -8,

examples/2D_tumbling_rectangle/case.py

@@ -60,6 +60,7 @@
             # Set IB to True and add 1 patch
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             "viscous": "T",
             # Formatted Database Files Structure Parameters
             "format": "silo",

examples/3D_ibm_bowshock/case.py

@@ -67,6 +67,7 @@
             # Set IB to True and add 1 patch
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             "viscous": "T",
             # Formatted Database Files Structure Parameters
             "format": "silo",

examples/3D_mibm_sphere_head_on_collision/case.py

@@ -83,6 +83,7 @@
             # Set IB to True and add 1 patch
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             # Formatted Database Files Structure Parameters
             "format": "silo",
             "precision": "double",

examples/3D_rotating_sphere/case.py

@@ -62,6 +62,7 @@
             # Set IB to True and add 1 patch
             "ib": "T",
             "num_ibs": 1,
+            "fd_order": 2,
             "viscous": "T",
             # Formatted Database Files Structure Parameters
             "format": "silo",

src/simulation/m_derived_variables.fpp

@@ -21,16 +21,7 @@ module m_derived_variables
     private; public :: s_initialize_derived_variables_module, s_initialize_derived_variables, s_compute_derived_variables, &
         & s_finalize_derived_variables_module
 
-    !> @name Finite-difference coefficients Finite-difference (fd) coefficients in x-, y- and z-coordinate directions. Note that
-    !! because sufficient boundary information is available for all the active coordinate directions, the centered family of the
-    !! finite-difference schemes is used.
-    !> @{
-    real(wp), public, allocatable, dimension(:,:) :: fd_coeff_x
-    real(wp), public, allocatable, dimension(:,:) :: fd_coeff_y
-    real(wp), public, allocatable, dimension(:,:) :: fd_coeff_z
-    !> @}
-
-    $:GPU_DECLARE(create='[fd_coeff_x, fd_coeff_y, fd_coeff_z]')
+    ! fd_coeff_x, fd_coeff_y, fd_coeff_z: declared in m_global_parameters so m_viscous can see them too
 
     ! @name Variables for computing acceleration
     !> @{
@@ -50,7 +41,7 @@ contains
         ! to be implemented in the subroutine s_compute_finite_difference_coefficients.
 
         ! Allocating centered finite-difference coefficients
-        if (probe_wrt) then
+        if (probe_wrt .or. ib) then
             @:ALLOCATE(fd_coeff_x(-fd_number:fd_number, 0:m))
             if (n > 0) then
                 @:ALLOCATE(fd_coeff_y(-fd_number:fd_number, 0:n))
@@ -74,9 +65,9 @@ contains
     !> Allocate and open derived variables. Computing FD coefficients.
     impure subroutine s_initialize_derived_variables
 
-        if (probe_wrt) then
+        if (probe_wrt .or. ib) then
             ! Opening and writing header of flow probe files
-            if (proc_rank == 0) then
+            if (proc_rank == 0 .and. probe_wrt) then
                 call s_open_probe_files()
                 call s_open_com_files()
             end if

src/simulation/m_global_parameters.fpp

@@ -173,6 +173,14 @@ module m_global_parameters
     integer :: fd_number  !< Finite-difference half-stencil size: MAX(1, fd_order/2)
     $:GPU_DECLARE(create='[fd_number]')
 
+    !> @name Centered finite-difference coefficients in x-, y- and z-coordinate directions
+    !> @{
+    real(wp), allocatable, dimension(:,:) :: fd_coeff_x
+    real(wp), allocatable, dimension(:,:) :: fd_coeff_y
+    real(wp), allocatable, dimension(:,:) :: fd_coeff_z
+    !> @}
+    $:GPU_DECLARE(create='[fd_coeff_x, fd_coeff_y, fd_coeff_z]')
+
     ! probe, integral: auto-generated in generated_decls.fpp
 
     !> @name Reference density and pressure for Tait EOS
@@ -818,15 +826,7 @@ contains
 
         if (ib) allocate (MPI_IO_IB_DATA%var%sf(0:m,0:n,0:p))
 
-        if (elasticity) then
-            fd_number = max(1, fd_order/2)
-        end if
-
-        if (mhd) then
-            fd_number = max(1, fd_order/2)
-        end if
-
-        if (probe_wrt) then
+        if (elasticity .or. mhd .or. probe_wrt .or. ib) then
             fd_number = max(1, fd_order/2)
         end if
 

src/simulation/m_ibm.fpp

@@ -910,9 +910,9 @@ contains
         integer                                                        :: i, j, k, l, encoded_ib_idx, ib_idx, ib_idx_temp, fluid_idx
         real(wp), dimension(num_ibs, 3)                                :: forces, torques
         ! viscous stress tensor with temp vectors to hold divergence calculations
-        real(wp), dimension(1:3,1:3) :: viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2
+        real(wp), dimension(1:3,1:3) :: viscous_stress
         real(wp), dimension(1:3)     :: local_force_contribution, radial_vector, local_torque_contribution
-        real(wp)                     :: cell_volume, dx, dy, dz, dynamic_viscosity
+        real(wp)                     :: cell_volume, dynamic_viscosity
 
         #:if not MFC_CASE_OPTIMIZATION and USING_AMD
             real(wp), dimension(3) :: dynamic_viscosities
@@ -937,8 +937,7 @@ contains
 
         $:GPU_PARALLEL_LOOP(private='[i, j, k, l, ib_idx, ib_idx_temp, encoded_ib_idx, fluid_idx, radial_vector, &
                             & local_force_contribution, cell_volume, local_torque_contribution, dynamic_viscosity, &
-                            & viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, dx, dy, dz]', copy='[forces, &
-                            & torques]', copyin='[dynamic_viscosities]', collapse=3)
+                            & viscous_stress]', copy='[forces, torques]', copyin='[dynamic_viscosities]', collapse=3)
         do i = 0, m
             do j = 0, n
                 do k = 0, p
@@ -948,34 +947,21 @@ contains
                         call s_get_neighborhood_idx(ib_idx_temp, ib_idx)  ! global patch ID -> local index
                         if (ib_idx > 0) then
                             ! get the vector pointing to the grid cell from the IB centroid
-                            if (num_dims == 3) then
-                                radial_vector = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, &
-                                                      & patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid]
-                            else
-                                radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, &
-                                                      & patch_ib(ib_idx)%y_centroid, 0._wp]
-                            end if
-                            dx = x_cc(i + 1) - x_cc(i)
-                            dy = y_cc(j + 1) - y_cc(j)
+                            radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, &
+                                                  & patch_ib(ib_idx)%y_centroid, 0._wp]
+                            if (num_dims == 3) radial_vector(3) = patch_ib(ib_idx)%z_centroid
 
                             local_force_contribution(:) = 0._wp
-                            do fluid_idx = 0, num_fluids - 1
-                                ! Get the pressure contribution to force via a finite difference to compute the 2D components of the
-                                ! gradient of the pressure and cell volume
-                                local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(eqn_idx%E &
-                                                         & + fluid_idx)%sf(i + 1, j, &
-                                                         & k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i - 1, j, k))/(2._wp*dx)  ! force is the negative pressure gradient
-                                local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(eqn_idx%E &
-                                                         & + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, &
-                                                         & j - 1, k))/(2._wp*dy)
-                                cell_volume = abs(dx*dy)
-                                ! add the 3D component of the pressure gradient, if we are working in 3 dimensions
+
+                            ! compute the pressure force component, which is the negative pressure gradient
+                            do l = -fd_number, fd_number
+                                local_force_contribution(1) = local_force_contribution(1) - (fd_coeff_x(l, &
+                                                         & i)*q_prim_vf(eqn_idx%E)%sf(i + l, j, k))
+                                local_force_contribution(2) = local_force_contribution(2) - (fd_coeff_y(l, &
+                                                         & j)*q_prim_vf(eqn_idx%E)%sf(i, j + l, k))
                                 if (num_dims == 3) then
-                                    dz = z_cc(k + 1) - z_cc(k)
-                                    local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(eqn_idx%E &
-                                                             & + fluid_idx)%sf(i, j, &
-                                                             & k + 1) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j, k - 1))/(2._wp*dz)
-                                    cell_volume = abs(cell_volume*dz)
+                                    local_force_contribution(3) = local_force_contribution(3) - (fd_coeff_z(l, &
+                                                             & k)*q_prim_vf(eqn_idx%E)%sf(i, j, k + l))
                                 end if
                             end do
 
@@ -989,41 +975,30 @@ contains
                                         & k)*dynamic_viscosities(fluid_idx))
                                 end do
 
-                                ! get the linear force components first
-                                call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, &
-                                                                     & j, k)
-                                call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, &
-                                                                     & j, k)
-                                ! get x derivative of the first-row of viscous stress tensor
-                                viscous_stress_div(1,1:3) = (viscous_stress_div_2(1,1:3) - viscous_stress_div_1(1,1:3))/(2._wp*dx)
-                                ! add the x components of the divergence to the force
-                                local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1,1:3)
-
-                                call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, &
-                                                                     & j - 1, k)
-                                call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, &
-                                                                     & j + 1, k)
-                                ! get y derivative of the second-row of viscous stress tensor
-                                viscous_stress_div(2,1:3) = (viscous_stress_div_2(2,1:3) - viscous_stress_div_1(2,1:3))/(2._wp*dy)
-                                ! add the y components of the divergence to the force
-                                local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2,1:3)
+                                do l = -fd_number, fd_number
+                                    call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i + l, j, k)
+                                    local_force_contribution(1:3) = local_force_contribution(1:3) + fd_coeff_x(l, &
+                                                             & i)*viscous_stress(1,1:3)
 
-                                if (num_dims == 3) then
-                                    call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, &
-                                                                         & j, k - 1)
-                                    call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, &
-                                                                         & j, k + 1)
-                                    viscous_stress_div(3,1:3) = (viscous_stress_div_2(3,1:3) - viscous_stress_div_1(3, &
-                                                       & 1:3))/(2._wp*dz)
-                                    ! add the z components of the divergence to the force
-                                    local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3,1:3)
-                                end if
+                                    call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i, j + l, k)
+                                    local_force_contribution(1:3) = local_force_contribution(1:3) + fd_coeff_y(l, &
+                                                             & j)*viscous_stress(2,1:3)
+
+                                    if (num_dims == 3) then
+                                        call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i, j, &
+                                                                             & k + l)
+                                        local_force_contribution(1:3) = local_force_contribution(1:3) + fd_coeff_z(l, &
+                                                                 & k)*viscous_stress(3,1:3)
+                                    end if
+                                end do
                             end if
 
                             call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution)
 
                             ! Update the force and torque values atomically to prevent race conditions
-                            do l = 1, 3
+                            cell_volume = dx(i)*dy(j)
+                            if (num_dims == 3) cell_volume = cell_volume*dz(k)
+                            do l = 1, num_dims
                                 $:GPU_ATOMIC(atomic='update')
                                 forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume)
                                 $:GPU_ATOMIC(atomic='update')