1104 mibm higher order central differencing#1614
Conversation
Claude Code ReviewHead SHA: a13521b Files changed:
Findings[Correctness] 3D radial-vector bug in 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_centroidThe initial array constructor correctly sets if (num_dims == 3) radial_vector(3) = z_cc(k) - patch_ib(ib_idx)%z_centroid[Memory]
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 ifBut if (probe_wrt) then
deallocate (accel_mag, x_accel)
...
end ifWhen |
|
Fixed the AI review comment |
…ell. I now get identical results to the python version
|
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 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. |
|
it certainly seems better, but remind me what we think accounts for the difference between the reference and our result |
|
@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. |
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1614 +/- ##
==========================================
- Coverage 60.49% 60.48% -0.01%
==========================================
Files 83 83
Lines 19902 19888 -14
Branches 2951 2951
==========================================
- Hits 12040 12030 -10
+ Misses 5868 5864 -4
Partials 1994 1994 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
@sbryngelson Did a higher resolution sim over night. This is 3000x3000 (100 grid cells across the cylinder diameter) of the Re=40 Ma=0.1 case MFC computed a drag coefficient of 1.549 compared to the reference value of 1.540 from 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 is worth documenting in an example and the readme since it's central to our claim that we have IBM working with the correct forces. |
|
otherwise it seems ready to merge |
|
the |
|
I don't think it is as clear-cut as that. When we run with Separating them gives more flexibility, but setting all 4 FD orders would instantiate a whole lot of memory that is somewhat wasted. The benefit of doing it the way I did was first that I did not add more case file parameters, but mostly that I am piggy-backing off already instantiated memory values. I think a better refactor would probably be drop the coefficient arrays altogether. They are extremely extremely cheap to compute, so storing them is wasteful. I feel like it is more efficient to compute the coefficient subroutine directly inside the loop to get the coefficient values. But doing this way would mean that we can have different finite difference orders for all 4 different cases in the same case file and that the net memory usage goes down. I do not know exactly how much it costs to compute them, but my intuition and basic reasoning tells me that it is almost zero. |
|
these input parameters are also used by post process in an entirely different manner, e.g., for schlieren i think. that's mostly what i mean. so we have a checker that checks if X postprocess parameter is true and if fd_order > 2 (for example). and setting fd_order the way you do conflicts with that. the main problem isn't your idea, it's that MFC already overloaded this parameter. re: compute at runtime vs. store in array. i think keeping in array might make sense because the array is small and used in many places. DRY. |
|
Opened the issue and copied the relevant comments. Discussion about this can move there. |
2 participants
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_forcessubroutineCloses #1104
Type of change (delete unused ones)
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.See the developer guide for full coding standards.
GPU changes (expand if you modified
src/simulation/)