1104 mibm higher order central differencing

benchmarks/ibm/case.py

@@ -68,6 +68,7 @@
             "riemann_solver": 2,
             "wave_speeds": 1,
             "avg_state": 2,
+            "fd_order": 2,
             "bc_x%beg": -3,
             "bc_x%end": -3,
             "bc_y%beg": -3,

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_ibm_viscous_drag_over_cylinder/case.py

@@ -0,0 +1,106 @@
+import json
+import math
+import sys
+
+"""
+A case file meant to perform a direct comparison with:
+Canuto D, Taira K. Two-dimensional compressible viscous flow around a circular cylinder. Journal of Fluid Mechanics. 2015;785:349-371. doi:10.1017/jfm.2015.635
+
+
+This performs a comparison to the low reynolds case in Table 4 of that paper
+"""
+
+Re_p = 40.0
+Ma_inf = 0.1
+d_cyl = 1.0
+r_cyl = d_cyl / 2.0
+
+gamma = 1.4
+U_inf = 1.0
+rho_inf = 1.0
+P_inf = rho_inf * U_inf**2 / (gamma * Ma_inf**2)
+
+case = {
+    # --- Output ---
+    "format": 1,
+    "precision": 2,
+    "run_time_info": "T",
+    "parallel_io": "T",
+    "prim_vars_wrt": "T",
+    "ib_state_wrt": "T",
+    # --- Domain ---
+    "x_domain%beg": -15.0,
+    "x_domain%end": 15.0,
+    "y_domain%beg": -15.0,
+    "y_domain%end": 15.0,
+    "m": 3000,
+    "n": 3000,
+    "p": 0,
+    "cyl_coord": "F",
+    # --- Time stepping ---
+    "cfl_adap_dt": "T",
+    "cfl_target": 0.5,
+    "n_start": 0,
+    "t_save": 10.0,
+    "t_stop": 100.0,
+    # --- Numerics ---
+    "num_patches": 1,
+    "num_fluids": 1,
+    "model_eqns": 2,
+    "alt_soundspeed": "F",
+    "mpp_lim": "F",
+    "mixture_err": "T",
+    "time_stepper": 3,
+    "weno_order": 5,
+    "weno_eps": 1.0e-10,
+    "weno_Re_flux": "T",
+    "weno_avg": "T",
+    "avg_state": 2,
+    "mapped_weno": "T",
+    "null_weights": "F",
+    "mp_weno": "F",
+    "riemann_solver": 2,
+    "low_Mach": 2,  # enable low mach disipation setting
+    "wave_speeds": 1,
+    "viscous": "T",
+    "fd_order": 4,
+    # --- Initial condition: uniform freestream, 2D rectangle ---
+    "patch_icpp(1)%geometry": 3,  # 2D rectangle
+    "patch_icpp(1)%x_centroid": 0.0,
+    "patch_icpp(1)%y_centroid": 0.0,
+    "patch_icpp(1)%length_x": 30.0,
+    "patch_icpp(1)%length_y": 30.0,
+    "patch_icpp(1)%vel(1)": U_inf,  # flow in +x
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%pres": P_inf,
+    "patch_icpp(1)%alpha_rho(1)": rho_inf,
+    "patch_icpp(1)%alpha(1)": 1.0,
+    # --- Fluid properties (calorically perfect air) ---
+    "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),  # MFC: 1/(gamma-1)
+    "fluid_pp(1)%pi_inf": 0.0,
+    "fluid_pp(1)%Re(1)": Re_p / (d_cyl * U_inf),
+    # --- Boundary conditions ---
+    "bc_x%beg": -7,
+    "bc_x%grcbc_in": "T",  # grbc for low-mach case
+    "bc_x%vel_in(1)": U_inf,
+    "bc_x%vel_in(2)": 0.0,
+    "bc_x%pres_in": P_inf,
+    "bc_x%alpha_rho_in(1)": rho_inf,
+    "bc_x%alpha_in(1)": 1.0,
+    "bc_x%end": -8,
+    "bc_x%grcbc_out": "T",
+    "bc_x%pres_out": P_inf,
+    "bc_y%beg": -9,
+    "bc_y%end": -9,
+    # --- IBM: 2D circle (cylinder cross-section, geometry=2) ---
+    "ib": "T",
+    "num_ibs": 1,
+    "patch_ib(1)%geometry": 2,
+    "patch_ib(1)%x_centroid": 0.0,
+    "patch_ib(1)%y_centroid": 0.0,
+    "patch_ib(1)%radius": r_cyl,
+    "patch_ib(1)%slip": "F",
+    "patch_ib(1)%moving_ibm": 0,
+}
+
+print(json.dumps(case, indent=4))

examples/2D_ibm_viscous_drag_over_cylinder/readme.md

@@ -0,0 +1,31 @@
+# 2D Viscous Drag Coefficient Convergence Test
+
+## Case Set Up
+
+This case file is an example of the convergence of the viscous drag over the surface of a 2D cylinder. We directly compare this result with the results of Canuto and Taira [1] which are the drag coefficient around a 2D cylinder in low-mach viscous flows for varying Mach and Reynolds numbers. The case set up is the exact parameters used by Taira and Colonius [2] with the additional determination of the freestream pressure from the freestream Mach number. A full list of coefficient and flow parameters can be determined in Table 4 of [1]. This case matches their setup for `Re = 40` and `Ma = 0.1`.
+
+## Numerics
+
+Of particular note in getting the numerics of this case to match closely to those of Canuto and Taira, we used two techniques:
+
+1. `"low_mach": 2` enables an improved WENO reconstruction scheme which reduces low-mach disipation from Thornbet et al. [3]
+
+2. `"mp_weno": "F"` disables MPWENO from Dinshaw, Balsara, and Shu [4] which clamps the WENO flux terms in an attempt to correct in the case of sharp discontinuities. However, these low-Mach low-Renolds cases are very smooth. So the clipping is not needed, but any triggerings could reduce our numeric accuracy. This proved to be a significant improvement.
+
+## Result
+
+To test the required convergence, comparisons to 4 different pressure and viscous drag force integration methods was made: surface integration and volume integration both computed at 2nd and 4th order. These values were compared in a post-processing step. In this case file, we make use of `fd_order: 4`, which enabled 4th order finite differencing techniques to be used in the IB drag force volume integraiton method of MFC. The following plot demonstrates that the MFC calculation matches the post-processing impliementation of 4th-order, and that this technique gives the appropriate convergence.
+
+![Screenshot](drag_coeffs.png)
+
+The final MFC-obtained value of the drag coefficient is $C_{d}^{MFC} = 1.549$, which differs only slightly from the value from Canuto and Taira of $C_{d}^{Canuto} = 1.540$.
+
+### References
+
+[1] Canuto D, Taira K. Two-dimensional compressible viscous flow around a circular cylinder. Journal of Fluid Mechanics. 2015;785:349-371. doi:10.1017/jfm.2015.635
+
+[2] Taira K, Colonius T. The immersed boundary method: A projection approach. Journal of Computational Physics. August 2007. doi:10.1016/j.jcp.2007.03.005
+
+[3] B. Thornber, A. Mosedale, D. Drikakis, D. Youngs, R.J.R. Williams, An improved reconstruction method for compressible flows with low Mach number features. Journal of Computational Physics. Volume 227, Issue 10,2008, Pages 4873-4894. doi:10.1016/j.jcp.2008.01.036.
+
+[4] Dinshaw S. Balsara, Chi-Wang Shu, Monotonicity Preserving Weighted Essentially Non-oscillatory Schemes with Increasingly High Order of Accuracy. Journal of Computational Physics. olume 160, Issue 2, 2000, Pages 405-452, doi:10.1006/jcph.2000.6443.

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