Fix IBM force: accumulate from both predictor AND corrected velocity

Cd was ~3x too low because only the predictor IBM call contributed to
force accumulation. At quasi-steady state u* inside body ≈ 0 (IBM zeroed
it last step), so predictor contribution is tiny. The missing force is
from the second IBM call: u^{n+1} = u*_IBM - dt*grad(p) inside body is
non-zero due to pressure correction; re-zeroing it gives the pressure drag.

Fix: add reset_force_accumulator() (called once per step in solver.cpp),
change both apply_forcing calls to ADD to last_Fx_/Fy_/Fz_. Total force
= predictor correction + pressure correction = full IBM drag.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: sbryngelson

include/ibm_forcing.hpp

@@ -69,6 +69,11 @@ class IBMForcing {
     /// Whether GPU data is mapped
     bool is_gpu_ready() const { return gpu_mapped_; }
 
+    /// Reset force accumulator to zero; call once at the start of each time step,
+    /// before the first apply_forcing call. Both IBM calls (predictor + corrected
+    /// velocity) then ADD their contributions so the total captures pressure drag too.
+    void reset_force_accumulator();
+
     /// Compute drag and lift forces on the body
     /// @param vel  Current velocity field
     /// @param dt   Time step size

src/ibm_forcing.cpp

@@ -254,6 +254,12 @@ void IBMForcing::unmap_from_gpu() {
     gpu_mapped_ = false;
 }
 
+void IBMForcing::reset_force_accumulator() {
+    last_Fx_ = 0.0;
+    last_Fy_ = 0.0;
+    last_Fz_ = 0.0;
+}
+
 void IBMForcing::apply_forcing(VectorField& vel, double dt) {
     const int Nx = mesh_->Nx;
     const int Ny = mesh_->Ny;
@@ -262,14 +268,13 @@ void IBMForcing::apply_forcing(VectorField& vel, double dt) {
     const bool is2D = mesh_->is2D();
     int Nz_eff = is2D ? 1 : Nz;
 
-    // Accumulate force on body from predictor velocity (before weight multiply)
-    // Only accumulate when dt > 0 (first call on velocity_star_); skip second call (dt=0).
+    // Accumulate IBM momentum correction into force accumulator (ADD, don't reset).
+    // Call reset_force_accumulator() once per step before the first apply_forcing call.
+    // Both calls (predictor + corrected velocity) contribute: predictor captures viscous
+    // damping, corrected velocity captures pressure drag (u^{n+1} = u*_IBM - dt*grad(p)).
     if (dt > 0.0) {
         const double dx = mesh_->dx;
         const double dz_val = is2D ? 1.0 : mesh_->dz;
-        last_Fx_ = 0.0;
-        last_Fy_ = 0.0;
-        last_Fz_ = 0.0;
 
         for (int k = Ng; k < Nz_eff + Ng; ++k) {
             for (int j = Ng; j < Ny + Ng; ++j) {
@@ -362,8 +367,9 @@ void IBMForcing::apply_forcing_device(double* u_ptr, double* v_ptr, double* w_pt
     const int v_n = static_cast<int>(v_total_);
     const int w_n = static_cast<int>(w_total_);
 
-    // Accumulate force on body from predictor velocity (before weight multiply)
-    // Only accumulate when dt > 0 (first call on velocity_star_); skip second call (dt=0).
+    // Accumulate IBM momentum correction (ADD to accumulator; caller resets via
+    // reset_force_accumulator() once per step). Both calls (predictor + corrected
+    // velocity) contribute so the total captures viscous + pressure drag.
     if (dt > 0.0) {
         [[maybe_unused]] const double dV = mesh_->dx * mesh_->dy * (mesh_->is2D() ? 1.0 : mesh_->dz);
         double Fx_acc = 0.0;
@@ -390,9 +396,9 @@ void IBMForcing::apply_forcing_device(double* u_ptr, double* v_ptr, double* w_pt
             }
         }
 
-        last_Fx_ = Fx_acc / dt * dV;
-        last_Fy_ = Fy_acc / dt * dV;
-        last_Fz_ = Fz_acc / dt * dV;
+        last_Fx_ += Fx_acc / dt * dV;
+        last_Fy_ += Fy_acc / dt * dV;
+        last_Fz_ += Fz_acc / dt * dV;
     }
 
     #pragma omp target teams distribute parallel for \

src/solver.cpp

@@ -1751,8 +1751,11 @@ double RANSSolver::step() {
     }
 
     // 3b. Apply IBM forcing to predicted velocity (before Poisson solve)
-    // Must apply to velocity_star_ (the predictor), not velocity_ (current step)
+    // Must apply to velocity_star_ (the predictor), not velocity_ (current step).
+    // Reset force accumulator here so both IBM calls (predictor + corrected velocity)
+    // contribute: predictor captures viscous damping, second call captures pressure drag.
     if (ibm_) {
+        ibm_->reset_force_accumulator();
         if (gpu_ready_ && ibm_->is_gpu_ready()) {
             ibm_->apply_forcing_device(velocity_star_u_ptr_, velocity_star_v_ptr_,
                                        mesh_->is2D() ? nullptr : velocity_star_w_ptr_,
@@ -2220,14 +2223,16 @@ double RANSSolver::step() {
     }
 
     // 5b. Re-apply IBM forcing after velocity correction
-    // Pressure correction may introduce non-zero velocity inside the body
-    // dt=0.0 (default): do not overwrite last_Fx_/Fy_/Fz_ accumulated from predictor
+    // Re-force corrected velocity: pressure correction may introduce non-zero velocity
+    // inside the body. Pass current_dt_ so this call also contributes to force accumulator
+    // (pressure drag u^{n+1} = u*_IBM - dt*grad(p) inside body captures pressure force).
     if (ibm_) {
         if (gpu_ready_ && ibm_->is_gpu_ready()) {
             ibm_->apply_forcing_device(velocity_u_ptr_, velocity_v_ptr_,
-                                       mesh_->is2D() ? nullptr : velocity_w_ptr_);
+                                       mesh_->is2D() ? nullptr : velocity_w_ptr_,
+                                       current_dt_);
         } else {
-            ibm_->apply_forcing(velocity_, 0.0);
+            ibm_->apply_forcing(velocity_, current_dt_);
         }
     }
 

tests/test_ibm_forcing.cpp

@@ -140,9 +140,10 @@ void test_force_computation() {
         }
     }
 
-    // apply_forcing with dt>0 accumulates forces from predictor velocity,
-    // then zeroes body cells. compute_forces returns the cached values.
+    // reset_force_accumulator() then apply_forcing(dt>0) accumulates forces.
+    // compute_forces returns the cached values.
     const double dt = 0.001;
+    ibm.reset_force_accumulator();
     ibm.apply_forcing(vel, dt);
     auto [Fx, Fy, Fz] = ibm.compute_forces(vel, dt);