update()

void fl::WaveSimulation2D_Real::update

(

)

Definition at line 257 of file wave_simulation_real.cpp.hpp.

                                   {
    i16 *curr = (whichGrid == 0 ? grid1.data() : grid2.data());
    i16 *next = (whichGrid == 0 ? grid2.data() : grid1.data());
 
    // Update horizontal boundaries.
    for (fl::size j = 0; j < height + 2; ++j) {
        if (mXCylindrical) {
            curr[j * stride + 0] = curr[j * stride + width];
            curr[j * stride + (width + 1)] = curr[j * stride + 1];
        } else {
            curr[j * stride + 0] = curr[j * stride + 1];
            curr[j * stride + (width + 1)] = curr[j * stride + width];
        }
    }
 
    // Update vertical boundaries.
    for (fl::size i = 0; i < width + 2; ++i) {
        curr[0 * stride + i] = curr[1 * stride + i];
        curr[(height + 1) * stride + i] = curr[height * stride + i];
    }
 
    i32 mCourantSq32 = static_cast<i32>(mCourantSq);
    // Hoist the lower-saturation bound — see WaveSimulation1D_Real::update()
    // for the rationale. Fold half-duplex into the per-cell clamp instead
    // of running a second pass over the grid (especially expensive on
    // PSRAM-backed grids).
    const i32 q15_min = mHalfDuplex ? 0 : -32768;
 
    // Update each inner cell. Per-row hoist: lift j*stride out of the inner
    // loop, set up row pointers to curr/next/above/below, and mark them
    // FL_RESTRICT_PARAM so the optimizer knows curr and next don't alias.
    // This gives the compiler a clean dependency picture across iterations
    // (it can vectorize / interleave loads without re-deriving the address).
    //
    // Outer-level branch on the stencil keeps the inner loop branchless;
    // the two kernels are otherwise identical (Q15 multiply, damping,
    // clamp). FivePoint is the backward-compatible default; the wrapper
    // class WaveSimulation2D auto-selects NinePointIsotropic at high
    // super-sample factors where the anisotropy of the 5-point stencil
    // becomes visually obvious.
    const bool nine_point = (mStencil == LaplacianStencil::NinePointIsotropic);
    for (fl::size j = 1; j <= height; ++j) {
        const fl::size row = j * stride;
        const i16* FL_RESTRICT_PARAM row_curr  = curr + row;
        const i16* FL_RESTRICT_PARAM row_above = curr + (row - stride);
        const i16* FL_RESTRICT_PARAM row_below = curr + (row + stride);
        i16*       FL_RESTRICT_PARAM row_next  = next + row;
        for (fl::size i = 1; i <= width; ++i) {
            const i32 c = row_curr[i];
            i32 laplacian;
            if (nine_point) {
                // 9-point isotropic Laplacian, scaled up by 6 to keep
                // everything in integer arithmetic:
                //   6 * lap = (NW+NE+SW+SE) + 4*(N+S+E+W) - 20*C
                // The /6 is folded into the term computation below so we
                // pay it once per cell rather than four times.
                const i32 diag = (i32)row_above[i - 1] + row_above[i + 1] +
                                 row_below[i - 1] + row_below[i + 1];
                const i32 nbr  = (i32)row_above[i] + row_below[i] +
                                 row_curr[i - 1] + row_curr[i + 1];
                laplacian = diag + (nbr << 2) - 20 * c;
            } else {
                // Standard 5-point Laplacian: N + S + E + W - 4*C.
                laplacian = (i32)row_curr[i + 1] + row_curr[i - 1] +
                            row_above[i] + row_below[i] - (c << 2);
            }
            // Promote to i64 before the multiply. With the 2D CFL clamp at
            // 0.5, the 5-point worst case product is ~4.3e9 and the
            // 9-point (scaled by 6, so |lap| ~6x larger) is ~2.6e10 — both
            // past i32 max (2.15e9). The i64 promote also acts as a safety
            // net if the clamp is ever regressed.
            i64 product = static_cast<i64>(mCourantSq32) * laplacian;
            i32 term;
            if (nine_point) {
                // Undo the x6 scaling on the 9-point Laplacian here.
                // Integer division by a compile-time-constant 6 is lowered
                // to a reciprocal-multiply on Cortex-M4+ and to a small
                // shift+add on AVR; either way it's well under the cost
                // of computing the Laplacian itself.
                term = static_cast<i32>((product >> 15) / 6);
            } else {
                term = static_cast<i32>(product >> 15);
            }
            // f = -next[index] + 2 * curr[index] + mCourantSq * laplacian.
            i32 f = -(i32)row_next[i] + (c << 1) + term;
 
            // Apply damping with the precomputed Q15 decay multiplier —
            // see the 1D update for the rationale.
            f = static_cast<i32>(
                (static_cast<i64>(f) * mDampDecayQ15) >> 15);
 
            // Clamp f into [q15_min, 32767] in a single step — subsumes
            // both the Q15 saturation clamp and the half-duplex zero pass.
            row_next[i] = static_cast<i16>(fl::clamp(f, q15_min, static_cast<i32>(32767)));
        }
    }
 
    // Swap the roles of the grids.
    whichGrid ^= 1;
}

References fl::clamp(), FL_RESTRICT_PARAM, grid1, grid2, height, mCourantSq, mDampDecayQ15, mHalfDuplex, mStencil, mXCylindrical, fl::NinePointIsotropic, stride, whichGrid, and width.

Here is the call graph for this function: