diff --git a/examples/2D_ibm_force_decomposition/README.md b/examples/2D_ibm_force_decomposition/README.md new file mode 100644 index 000000000..b133f5b33 --- /dev/null +++ b/examples/2D_ibm_force_decomposition/README.md @@ -0,0 +1,43 @@ +# Immersed-boundary force under a changing MPI decomposition + +A static cylinder in uniform flow at Re 40, Ma 0.1, on a uniform 200 x 200 grid. The body sits at the +origin, which is also the center of the domain, so any decomposition with an even number of ranks in a +direction puts a subdomain edge straight through it. + +The force on a rigid body is a property of the flow. Running the same case on a different number of ranks +must not change it. This case measures whether it does. + +## Running it + +``` +./mfc.sh run examples/2D_ibm_force_decomposition/case.py -n 1 +./mfc.sh run examples/2D_ibm_force_decomposition/case.py -n 4 +``` + +Each run writes `restart_data/ib_state_200.dat`: 20 doubles, `[time, Fx, Fy, Fz, Tx, Ty, Tz, ...]`. + +``` +python3 -c "import numpy as np; a = np.fromfile('restart_data/ib_state_200.dat'); print(a[1], a[2])" +``` + +## What it shows + +`fd_order = 4` here, so the force integral's stencil reaches two cells, and a body cell two cells from a +subdomain edge asks for finite-difference coefficients outside the interior. Before the coefficients were +defined over that range they were read off the end of the array: + +| ranks | Fx before | Fx after | +| --- | --- | --- | +| 1 | 1.02352019 | 1.02352019 | +| 4 | 1.02933438 | 1.02351629 | + +0.57 percent of drag appearing out of adjacent memory, against 0.0004 percent after. The single-rank +answer is unchanged, because a single rank never reaches outside its own interior. + +The size of the discrepancy is set by whatever is adjacent in memory and is not bounded by anything: on a +3D sphere at Re 100 on 64 ranks the same read produced a transverse force of 1.08 times the drag on a body +that has none. + +Lift is a second, independent check: the cylinder is symmetric about y = 0, so Fy must be zero. It is +about 8e-6 here in every configuration, which is the level at which the staircase representation of a +circle on this grid is symmetric, and it does not move. diff --git a/examples/2D_ibm_force_decomposition/case.py b/examples/2D_ibm_force_decomposition/case.py new file mode 100644 index 000000000..26865fb4a --- /dev/null +++ b/examples/2D_ibm_force_decomposition/case.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +"""The smallest case that shows MFC's immersed-boundary force depending on the MPI decomposition. + +A static cylinder in uniform flow, Re 40, on a grid coarse enough to run on a laptop. The body sits at the +origin, which is also the domain center, so a decomposition with an even number of ranks in a direction puts +a subdomain edge straight through the body. Run it at one rank and at four and compare `restart_data/ib_state_200.dat`: +the cylinder is symmetric about y = 0 and the lift must be zero, and the drag must not care how the domain +was cut up. +""" + +import json +import os + +Re, Ma, d, gamma, U, rho = 40.0, 0.1, 1.0, 1.4, 1.0, 1.0 +P = rho * U**2 / (gamma * Ma**2) +L = float(os.environ.get("L", 5.0)) +dx = float(os.environ.get("DX", 0.1)) +N = int(2 * L / dx) +NSTEP = int(os.environ.get("NSTEP", 200)) + +case = { + "run_time_info": "T", + "parallel_io": "T", + "prim_vars_wrt": "T", + "ib_state_wrt": "T", + "format": "silo", + "precision": "double", + "x_domain%beg": -L, + "x_domain%end": L, + "y_domain%beg": -L, + "y_domain%end": L, + "m": N - 1, + "n": N - 1, + "p": 0, + "cyl_coord": "F", + "dt": 0.3 * dx / (U + U / Ma), + "t_step_start": 0, + "t_step_stop": NSTEP, + "t_step_save": NSTEP, + "num_patches": 1, + "num_fluids": 1, + "model_eqns": "5eq", + "alt_soundspeed": "F", + "mpp_lim": "F", + "mixture_err": "T", + "time_stepper": "rk3", + "weno_order": 5, + "weno_eps": 1.0e-10, + "weno_Re_flux": "T", + "weno_avg": "T", + "avg_state": "arithmetic", + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": "hllc", + "low_Mach": 2, + "wave_speeds": "direct", + "viscous": "T", + "fd_order": int(os.environ.get("FDORDER", 4)), + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 0.0, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%length_x": 2 * L, + "patch_icpp(1)%length_y": 2 * L, + "patch_icpp(1)%vel(1)": U, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%pres": P, + "patch_icpp(1)%alpha_rho(1)": rho, + "patch_icpp(1)%alpha(1)": 1.0, + "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0), + "fluid_pp(1)%eos": "ideal_gas", + "fluid_pp(1)%Re(1)": Re / (d * U), + "bc_x%beg": -7, + "bc_x%grcbc_in": "T", + "bc_x%vel_in(1)": U, + "bc_x%vel_in(2)": 0.0, + "bc_x%pres_in": P, + "bc_x%alpha_rho_in(1)": rho, + "bc_x%alpha_in(1)": 1.0, + "bc_x%end": -8, + "bc_y%beg": -8, + "bc_y%end": -8, + "ib": "T", + "num_ibs": 1, + "ib_neighborhood_radius": 3, + "patch_ib(1)%geometry": 2, + "patch_ib(1)%x_centroid": 0.0, + "patch_ib(1)%y_centroid": 0.0, + "patch_ib(1)%radius": d / 2, + "patch_ib(1)%slip": "F", + "patch_ib(1)%moving_ibm": 0, +} +print(json.dumps(case, indent=4)) diff --git a/src/simulation/m_derived_variables.fpp b/src/simulation/m_derived_variables.fpp index 6a4b92ddc..3569b2fdc 100644 --- a/src/simulation/m_derived_variables.fpp +++ b/src/simulation/m_derived_variables.fpp @@ -38,14 +38,16 @@ contains ! higher than fourth-order accuracy coefficients are wanted, the formulae required to compute these coefficients will have ! to be implemented in the subroutine s_compute_finite_difference_coefficients. - ! Allocating centered finite-difference coefficients + ! Allocating centered finite-difference coefficients. The IB force integral evaluates the viscous stress + ! at every cell its own stencil reaches, so it asks for coefficients up to fd_number outside the interior + ! whenever a body touches a subdomain edge; that range is included here rather than read off the end. if (probe_wrt .or. ib) then - @:ALLOCATE(fd_coeff_x(-fd_number:fd_number, 0:m)) + @:ALLOCATE(fd_coeff_x(-fd_number:fd_number, -fd_number:m + fd_number)) if (n > 0) then - @:ALLOCATE(fd_coeff_y(-fd_number:fd_number, 0:n)) + @:ALLOCATE(fd_coeff_y(-fd_number:fd_number, -fd_number:n + fd_number)) end if if (p > 0) then - @:ALLOCATE(fd_coeff_z(-fd_number:fd_number, 0:p)) + @:ALLOCATE(fd_coeff_z(-fd_number:fd_number, -fd_number:p + fd_number)) end if @:ALLOCATE(accel_mag(0:m, 0:n, 0:p)) @@ -63,6 +65,8 @@ contains !> Allocate and open derived variables. Computing FD coefficients. impure subroutine s_initialize_derived_variables + type(int_bounds_info) :: fd_offset + if (probe_wrt .or. ib) then ! Opening and writing header of flow probe files if (proc_rank == 0 .and. probe_wrt) then @@ -70,15 +74,16 @@ contains call s_open_com_files() end if ! Computing centered finite difference coefficients - call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x, buff_size, fd_number, fd_order) + fd_offset%beg = fd_number; fd_offset%end = fd_number + call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x, buff_size, fd_number, fd_order, fd_offset) $:GPU_UPDATE(device='[fd_coeff_x]') if (n > 0) then - call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y, buff_size, fd_number, fd_order) + call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y, buff_size, fd_number, fd_order, fd_offset) $:GPU_UPDATE(device='[fd_coeff_y]') end if if (p > 0) then - call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z, buff_size, fd_number, fd_order) + call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z, buff_size, fd_number, fd_order, fd_offset) $:GPU_UPDATE(device='[fd_coeff_z]') end if end if