diff --git a/examples/2D_backward_facing_step/case.py b/examples/2D_backward_facing_step/case.py index a06d9b7ac6..e7163d48d0 100644 --- a/examples/2D_backward_facing_step/case.py +++ b/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", diff --git a/examples/2D_forward_facing_step/case.py b/examples/2D_forward_facing_step/case.py index b7733a923e..c9bbde1289 100644 --- a/examples/2D_forward_facing_step/case.py +++ b/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", diff --git a/examples/2D_ibm/case.py b/examples/2D_ibm/case.py index a0d4eccb97..ea1799f946 100644 --- a/examples/2D_ibm/case.py +++ b/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", diff --git a/examples/2D_ibm_airfoil/case.py b/examples/2D_ibm_airfoil/case.py index 8afd9883ed..3aa7ee1588 100644 --- a/examples/2D_ibm_airfoil/case.py +++ b/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", diff --git a/examples/2D_ibm_cfl_dt/case.py b/examples/2D_ibm_cfl_dt/case.py index 434a527e8f..778af5d6a1 100644 --- a/examples/2D_ibm_cfl_dt/case.py +++ b/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", diff --git a/examples/2D_ibm_ellipse/case.py b/examples/2D_ibm_ellipse/case.py index 357f122319..cc82c3a275 100644 --- a/examples/2D_ibm_ellipse/case.py +++ b/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", diff --git a/examples/2D_ibm_multiphase/case.py b/examples/2D_ibm_multiphase/case.py index f4239b4f97..46a3749823 100644 --- a/examples/2D_ibm_multiphase/case.py +++ b/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", diff --git a/examples/2D_ibm_stl_MFCCharacter/case.py b/examples/2D_ibm_stl_MFCCharacter/case.py index 9616b97ef6..4a217a3840 100644 --- a/examples/2D_ibm_stl_MFCCharacter/case.py +++ b/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 diff --git a/examples/2D_ibm_stl_test/case.py b/examples/2D_ibm_stl_test/case.py index 843c6a535d..52efc63486 100644 --- a/examples/2D_ibm_stl_test/case.py +++ b/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 diff --git a/examples/2D_ibm_stl_wedge/case.py b/examples/2D_ibm_stl_wedge/case.py index a644ff1e0e..802502c7e6 100644 --- a/examples/2D_ibm_stl_wedge/case.py +++ b/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 diff --git a/examples/2D_mibm_cylinder_in_cross_flow/case.py b/examples/2D_mibm_cylinder_in_cross_flow/case.py index 0c89ef5962..bd05c60841 100644 --- a/examples/2D_mibm_cylinder_in_cross_flow/case.py +++ b/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, diff --git a/examples/2D_mibm_particle_cloud/case.py b/examples/2D_mibm_particle_cloud/case.py index 97c14b81d0..c6aafb6eb8 100644 --- a/examples/2D_mibm_particle_cloud/case.py +++ b/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, diff --git a/examples/2D_mibm_shock_cylinder/case.py b/examples/2D_mibm_shock_cylinder/case.py index f344ef9e27..b615e8afd4 100644 --- a/examples/2D_mibm_shock_cylinder/case.py +++ b/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, diff --git a/examples/2D_tumbling_rectangle/case.py b/examples/2D_tumbling_rectangle/case.py index ac7f16158f..5237f21258 100644 --- a/examples/2D_tumbling_rectangle/case.py +++ b/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", diff --git a/examples/3D_ibm_bowshock/case.py b/examples/3D_ibm_bowshock/case.py index 05c04d3917..ae3ce4b5d9 100644 --- a/examples/3D_ibm_bowshock/case.py +++ b/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", diff --git a/examples/3D_mibm_sphere_head_on_collision/case.py b/examples/3D_mibm_sphere_head_on_collision/case.py index 981aacb310..ee5bdbb2b4 100644 --- a/examples/3D_mibm_sphere_head_on_collision/case.py +++ b/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", diff --git a/examples/3D_rotating_sphere/case.py b/examples/3D_rotating_sphere/case.py index 1093211f59..9b1e5373d3 100644 --- a/examples/3D_rotating_sphere/case.py +++ b/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", diff --git a/src/simulation/m_derived_variables.fpp b/src/simulation/m_derived_variables.fpp index 22c0477065..6fc1518b58 100644 --- a/src/simulation/m_derived_variables.fpp +++ b/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 diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 8378324661..a51ffd897b 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/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 diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 811db7e54f..24507f0b6c 100644 --- a/src/simulation/m_ibm.fpp +++ b/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') diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp index d208428900..65d78db841 100644 --- a/src/simulation/m_viscous.fpp +++ b/src/simulation/m_viscous.fpp @@ -1263,34 +1263,28 @@ contains real(wp), intent(in) :: dynamic_viscosity integer, intent(in) :: i, j, k real(wp), dimension(1:3,1:3) :: velocity_gradient_tensor - real(wp), dimension(1:3) :: dx real(wp) :: divergence real(wp) :: mu_eff, gamma_dot_c integer :: l, q !< iterators integer :: fl + integer :: r - ! zero the viscous stress, collection of velocity derivatives, and spatial finite differences + ! zero the viscous stress and collection of velocity derivatives viscous_stress_tensor = 0._wp velocity_gradient_tensor = 0._wp - dx = 0._wp - ! get the change in x used in the finite difference equation - dx(1) = 0.5_wp*(x_cc(i + 1) - x_cc(i - 1)) - dx(2) = 0.5_wp*(y_cc(j + 1) - y_cc(j - 1)) - if (num_dims == 3) then - dx(3) = 0.5_wp*(z_cc(k + 1) - z_cc(k - 1)) - end if - - ! compute the velocity gradient tensor + ! compute the velocity gradient tensor with the same fd_order-respecting stencil as the stress-divergence outer derivative do l = 1, num_dims - velocity_gradient_tensor(l, 1) = (q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i + 1, j, & - & k) - q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i - 1, j, k))/(2._wp*dx(1)) - velocity_gradient_tensor(l, 2) = (q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j + 1, & - & k) - q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j - 1, k))/(2._wp*dx(2)) - if (num_dims == 3) then - velocity_gradient_tensor(l, 3) = (q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j, & - & k + 1) - q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j, k - 1))/(2._wp*dx(3)) - end if + do r = -fd_number, fd_number + velocity_gradient_tensor(l, 1) = velocity_gradient_tensor(l, 1) + fd_coeff_x(r, & + & i)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i + r, j, k) + velocity_gradient_tensor(l, 2) = velocity_gradient_tensor(l, 2) + fd_coeff_y(r, & + & j)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j + r, k) + if (num_dims == 3) then + velocity_gradient_tensor(l, 3) = velocity_gradient_tensor(l, 3) + fd_coeff_z(r, & + & k)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j, k + r) + end if + end do end do ! Non-Newtonian: per-sample mixture viscosity from the local strain rate, so each diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 5f49491883..ca631ba0b6 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -605,6 +605,8 @@ def check_ibm(self): ib_state_wrt = self.get("ib_state_wrt", "F") == "T" + fd_order = self.get("fd_order") + self.prohibit(ib and fd_order is None, "fd_order must be specified for ib") self.prohibit(ib and n <= 0, "Immersed Boundaries do not work in 1D (requires n > 0)") has_particle_clouds = num_particle_clouds > 0 and any((self.get(f"particle_cloud({i})%num_particles", 0) or 0) > 0 for i in range(1, num_particle_clouds + 1)) self.prohibit( diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index 495ae19d8b..a981e88edc 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -679,6 +679,7 @@ def alter_num_fluids(dimInfo): "fluid_pp(1)%nn": 0.5, "ib": "T", "num_ibs": 1, + "fd_order": 2, "ib_state_wrt": "T", "patch_ib(1)%geometry": 3, "patch_ib(1)%x_centroid": 0.5, @@ -876,6 +877,7 @@ def alter_ppn(dimInfo): "p": 49, "ib": "T", "num_ibs": 1, + "fd_order": 2, "patch_ib(1)%geometry": 8, "patch_ib(1)%x_centroid": 0.5, "patch_ib(1)%y_centroid": 0.5, @@ -901,6 +903,7 @@ def alter_ib(dimInfo, six_eqn_model=False, viscous=False): { "ib": "T", "num_ibs": 1, + "fd_order": 2, "patch_ib(1)%x_centroid": 0.5, "patch_ib(1)%y_centroid": 0.5, "patch_ib(1)%radius": 0.1, @@ -987,6 +990,7 @@ def alter_ib(dimInfo, six_eqn_model=False, viscous=False): { "ib": "T", "num_ibs": 1, + "fd_order": 2, "bc_x%beg": -1, "bc_x%end": -1, "bc_y%beg": -1, @@ -1008,6 +1012,7 @@ def ibm_stl(): common_mods = { "t_step_stop": Nt, "t_step_save": Nt, + "fd_order": 2, "num_stl_models": 1, "patch_ib(1)%model_id": 1, "stl_models(1)%model_scale(1)": 5.0,