1104 mibm higher order central differencing

[danieljvickers]

Description

For a while I have been putting off adding higher-order central difference mostly because I did not see an immediate need. However, it appears that the current viscous force calculation in MFC has some significant errors for low Reynolds cases. Increasing the finite-difference order for the derivatives to 4th order appears to significantly reduce this error back within acceptable tolerances. This branch integrates the IB code with the existing finite-difference coefficients in MFC to be utilized in the force calculation loop. This has added the side-benefit of somewhat reducing the amount of code in the s_compute_ib_forces subroutine

Closes #1104

Type of change (delete unused ones)

• Bug fix
• New feature
• Refactor

Testing

I did a significant amount of post-processing analysis of the data to determine the best integration method for low-Reynolds cases, which wound up being

Checklist

Check these like this [x] to indicate which of the below applies.

• I added or updated tests for new behavior
• I updated documentation if user-facing behavior changed

See the developer guide for full coding standards.

GPU changes (expand if you modified src/simulation/)

• GPU results match CPU results
• Tested on NVIDIA GPU or AMD GPU

[github-actions]

Claude Code Review

Head SHA: a13521b

Files changed:

• 7
• examples/2D_mibm_cylinder_in_cross_flow/case.py
• examples/2D_mibm_particle_cloud/case.py
• examples/2D_mibm_shock_cylinder/case.py
• src/simulation/m_derived_variables.fpp
• src/simulation/m_global_parameters.fpp
• src/simulation/m_ibm.fpp
• src/simulation/m_viscous.fpp

Findings

[Correctness] 3D radial-vector bug in m_ibm.fpp

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

The initial array constructor correctly sets radial_vector(1:2) to the cell-to-centroid offset. However, in 3D the last line replaces radial_vector(3) with z_centroid alone, instead of z_cc(k) - patch_ib(ib_idx)%z_centroid. The vector from the IB centroid to the grid cell is (x_cc(i) - x_centroid, y_cc(j) - y_centroid, z_cc(k) - z_centroid); the new code gives (x_cc(i) - x_centroid, y_cc(j) - y_centroid, z_centroid). This produces wrong torques for all 3D IB cases. The correct replacement is:

if (num_dims == 3) radial_vector(3) = z_cc(k) - patch_ib(ib_idx)%z_centroid

---

[Memory] accel_mag/x_accel/y_accel/z_accel allocated under ib but freed only under probe_wrt in m_derived_variables.fpp

s_initialize_derived_variables_module now allocates those four arrays whenever probe_wrt .or. ib:

if (probe_wrt .or. ib) then
    ...
    @:ALLOCATE(accel_mag(0:m, 0:n, 0:p))
    @:ALLOCATE(x_accel(0:m, 0:n, 0:p))
    ...
end if

But s_finalize_derived_variables_module (unchanged) frees them only when probe_wrt:

if (probe_wrt) then
    deallocate (accel_mag, x_accel)
    ...
end if

When ib .and. .not. probe_wrt these arrays are allocated, never used (the only consumer is inside if (probe_wrt) in s_compute_derived_variables), and never freed. The finalize condition should be widened to probe_wrt .or. ib, or the accel allocations should remain guarded by probe_wrt alone and be split out of the shared if block.

[danieljvickers]

Fixed the AI review comment

[danieljvickers]

I reran a convergence test and compared the force integration to my independent python implementation of all of the drag coefficient for on a 900x900 grid with Re=40 Mach 0.5 cylinder using fd_order: 4:

You can see that we now obtain the 4th-order integral. I take this as now working, and will proceed to updating tests before calling this complete.

[sbryngelson]

it certainly seems better, but remind me what we think accounts for the difference between the reference and our result

[danieljvickers]

@sbryngelson This is only on a 900x900 grid. I get much less error as I converge towards larger grid sizes in this case. This is to demonstrate that it matches the python solver that I have been using for analysis. I have separate results in my slides showing the accuracy of the python solver. Rerunning at a much higher resolution will just take a few hours, but I can do that if you want me to include it before merge.