Wave_Equation_Improvement

2025-11-19 19:41 UTC gpt-5-1 Open in ChatGPT ↗

Here’s what I can see and what I’d fix/improve based on the compute shader you shared.

1. Overall structure

You’ve got:

#pragma kernel CSMain
#pragma kernel Combine

Texture2D<float> PState;
Texture2D<float> State;
RWTexture2D<float> NState;

Texture2D<float4> Sources;
Texture2D<float> Obstacles;

float dt;
float c;
float D;
float _Time;

Conceptually this is a classic wave solver:

PState – previous timestep
State – current timestep
NState – next timestep (RW)
Sources – excitation / emitters
Obstacles – walls/geometry
dt, c, D – timestep, wave speed, damping
CSMain – presumably does the stencil / wave equation update
Combine – adds sources on top of the propagated state

That’s all reasonable. Most of the issues are in the details: types, indexing, and a couple of Unity/HLSL gotchas.

2. Biggest correctness issues

2.1. Reading from `Texture2D` with `[]`

In your Combine kernel:

[numthreads(8,8,1)]
void Combine (uint3 id : SV_DispatchThreadID)
{
    NState[id.xy] = State[id.xy] + Sources[id.xy];
}

Problem: State and Sources are Texture2D<...>, which are SRVs. In HLSL:

RWTexture2D<T> can be indexed with [] using integer coordinates.
Texture2D<T> cannot; you must call .Load() or .Sample().

So State[id.xy] and Sources[id.xy] are invalid and will not compile (or will behave incorrectly depending on compiler).

Fix options:

Option A – make them RW textures and use []:

RWTexture2D<float> PState;
RWTexture2D<float> State;
RWTexture2D<float> NState;
RWTexture2D<float4> Sources;
RWTexture2D<float> Obstacles;

Then State[id.xy] and Sources[id.xy] are fine.

Option B – keep them as Texture2D and use Load:

float stateVal   = State.Load(int3(id.xy, 0));
float4 srcVal    = Sources.Load(int3(id.xy, 0));
NState[id.xy]    = stateVal + srcVal.r; // or some combination of (r,g,b,a)

This is usually what you want if you are binding them as SRVs from the CPU side.

2.2. `Sources` is `float4` but `NState` is `float`

Right now:

Texture2D<float4> Sources;
RWTexture2D<float> NState;

NState[id.xy] = State[id.xy] + Sources[id.xy];

Assuming you fix the indexing, you still have a type mismatch:

State is float
Sources is float4
float + float4 → float4 (scalar is promoted to vector)
Assigning float4 to RWTexture2D<float> is not valid.

You must decide how you want to use Sources:

If it’s a greyscale mask: make it Texture2D<float> and keep everything scalar.
If you only care about one channel: explicitly pick that channel.

Suggested fix:

Texture2D<float> Sources; // Greyscale, simplest
// or:
float src = Sources.Load(int3(id.xy, 0)).r;

Then:

float stateVal = State.Load(int3(id.xy, 0));
float srcVal   = Sources.Load(int3(id.xy, 0));
NState[id.xy]  = stateVal + srcVal;

2.3. `_Time` naming conflict

Unity already defines a _Time uniform (typically float4 _Time).
You’re declaring:

float _Time;

This can cause:

name collisions
confusing behavior depending on include order and compiler

Fix: rename your uniform to something else, e.g.:

float SimulationTime;

Or, if you really want Unity’s _Time, do not redeclare it and instead use the built-in.

2.4. Likely missing bounds check

If you dispatch your compute shader with dimensions rounded up (common pattern in Unity: ceil(width / groupSize)), SV_DispatchThreadID may extend beyond the valid texture region.

Your kernels should guard against that:

cbuffer SimulationParams
{
    uint Width;
    uint Height;
    // dt, c, D, etc...
};

[numthreads(8,8,1)]
void Combine (uint3 id : SV_DispatchThreadID)
{
    if (id.x >= Width || id.y >= Height) return;

    float stateVal = State.Load(int3(id.xy, 0));
    float srcVal   = Sources.Load(int3(id.xy, 0));
    NState[id.xy]  = stateVal + srcVal;
}

Otherwise you risk out-of-bounds reads/writes if your dispatch is not perfectly aligned.

3. Wave equation / physics considerations

I can’t see your CSMain body (it’s redacted/omitted as ...), but for a typical 2D wave equation with damping you’d expect something like:

// Laplacian using 4-neighbour or 8-neighbour stencil
float lap =
    State.Load(int3(id.x + 1, id.y, 0)) +
    State.Load(int3(id.x - 1, id.y, 0)) +
    State.Load(int3(id.x, id.y + 1, 0)) +
    State.Load(int3(id.x, id.y - 1, 0)) -
    4.0 * State.Load(int3(id.xy, 0));

// Wave update (simple form)
float uPrev = PState.Load(int3(id.xy, 0));
float u     = State.Load(int3(id.xy, 0));

float dt2c2 = c * c * dt * dt;

float next = 2.0 * u - uPrev + dt2c2 * lap;

// Damping (e.g. simple linear term)
next *= exp(-D * dt);
NState[id.xy] = next;

Points to watch:

Stability: For an undamped wave equation on a regular grid, c * dt / dx must satisfy a CFL condition (e.g. ≤ 1/√2 for 2D with central differences).
If dt and c are too large, the simulation blows up.

Obstacles: Obstacles is declared but not shown in use. Usually you’d:

Zero out next where obstacles exist, or

Treat them as fixed boundary conditions, e.g.:

float obs = Obstacles.Load(int3(id.xy, 0));
if (obs > 0.5)
{
    next = 0.0;       // rigid wall
    // optionally: set PState = State = 0 here via another pass
}

Boundary conditions at edges:
- You need either clamped indexing, mirrored boundaries, or explicit conditions (e.g. zero / reflective).
- Make sure your Laplacian doesn’t read outside the texture; handle x=0, x=Width-1, etc.

4. API / Unity-side binding consistency

Double-check that your C# code matches the shader:

Resource types:
- If you’re binding with SetTexture(kernel, "State", stateRT), that RenderTexture must have enableRandomWrite set true for RWTexture2D, but not required for plain Texture2D SRVs.
- Decide which of your textures are read-only vs read-write, and ensure both the C# side and HLSL declarations line up.

Uniforms:

Use either a cbuffer or loose uniforms; for compute shaders, both are acceptable but cbuffer is often clearer:

cbuffer WaveParams
{
    float dt;
    float c;
    float D;
    float SimulationTime;
    uint Width;
    uint Height;
};

Thread group / dispatch:
- With [numthreads(8,8,1)], choose dispatch counts as:
```
int groupsX = (width  + 7) / 8;
int groupsY = (height + 7) / 8;
cs.Dispatch(kernel, groupsX, groupsY, 1);
```
- That matches the bounds check pattern mentioned earlier.

5. Small style / readability suggestions

Group configuration data:

Put all scalar params into a cbuffer for clarity and better packing.
Consistent naming:
- PState, State, NState is fine; sometimes Prev, Curr, Next is a bit clearer.
- Avoid leading underscores in custom variables in Unity compute shaders to reduce confusion with built-ins.
Comment data ranges:

Note in comments whether fields are normalized (e.g. height in [-1,1] vs [0,1]), what D physically represents (damping coefficient), etc. Helps future you tune them.

Quick summary of concrete fixes

Fix texture access:
- Use .Load(int3(id.xy, 0)) for Texture2D or convert to RWTexture2D if you want [] indexing.
Fix Sources type:
- Either use Texture2D<float> or explicitly use Sources.Load(...).r etc.
Rename _Time:
- Avoid clashing with Unity’s built-in _Time.
Add bounds checks and size uniforms if you dispatch with padded group counts.
Verify wave update uses a stable scheme (CFL condition) and that Obstacles + edges are handled.

If you paste the full CSMain body (the part that’s currently ... in the file I see), I can walk through the numerical scheme step-by-step and check stability, artifacts, and optimization opportunities (e.g. using groupshared for neighbor fetches).

Yeah, I know exactly what you mean 🙂

Shallow-water waves have a phase speed that depends on depth:

$ c = \sqrt{g,h} $

$h$ = water depth
$g$ = gravity

So:

Far from shore (deep water) → $h$ large → $c$ large → waves move faster.
Near shore (shallow) → $h$ small → $c$ small → waves slow down.

If a wavefront hits the coastline at an angle, one “side” of the wave will get into shallow water first and slow down first. The deeper side keeps going faster a bit longer. That difference in speed causes the wavefront to pivot/rotate, bending so that it becomes more and more parallel to the depth contours (and usually roughly parallel to the shoreline).

It’s literally refraction in the Snell’s-law sense, just with $c = \sqrt{g h}$ instead of light speeds in two media:

$ \frac{\sin\theta_1}{c_1} = \frac{\sin\theta_2}{c_2} $

Same idea: part of the front in “slower” medium bends toward the normal.

If you want, next step we can make your compute shader do this by:

Giving it a depth map.
Making your wave speed $c(x,y) \propto \sqrt{h(x,y)}$.
Letting the varying $c$ naturally bend the waves as they propagate (or even biasing propagation along depth gradients for extra control).

Nice, this is exactly the kind of fun tweak that makes the sim feel “smart” instead of just noisy 🙂

Let’s wire your solver up to a depth map–driven wave speed so refraction falls out naturally.

1. Add a depth map + spatially varying wave speed

We’ll:

Add a Depth texture.
Map depth → local wave speed $c(x,y)$.
Use $c(x,y)^2$ in the update step instead of a single global c.

I’ll show a full, drop-in version of your compute shader with those changes, staying close to your original naming.

Updated resource + parameter declarations

// Each #kernel tells which function to compile; you can have many kernels
#pragma kernel CSMain
#pragma kernel Combine

// States
Texture2D<float>   PState;     // previous state
Texture2D<float>   State;      // current state
RWTexture2D<float> NState;     // next state

// Inputs
Texture2D<float>   Sources;    // changed from float4 -> float
Texture2D<float>   Obstacles;  // 0 = water, 1 = solid
Texture2D<float>   Depth;      // depth map (0..1 or actual depth)

// Simulation constants
cbuffer WaveParams
{
    float dt;          // timestep
    float cShallow;    // wave speed in shallowest areas
    float cDeep;       // wave speed in deepest areas
    float D;           // damping
    uint  Width;       // simulation width
    uint  Height;      // simulation height
};

Depth → speed mapping:
We’ll treat Depth in [0,1], where 0 is shallow, 1 deep, and interpolate c(x,y) between cShallow and cDeep.

float GetLocalSpeed(int3 coord)
{
    float d = Depth.Load(coord);              // 0..1
    // map 0..1 depth to speed (shallow slower, deep faster)
    return lerp(cShallow, cDeep, d);
}

You can set cShallow and cDeep from C# to tune how dramatic the refraction is.

2. Depth-aware wave propagation (CSMain)

Here’s a complete CSMain using a 4-neighbour Laplacian and your depth-controlled c(x,y):

[numthreads(8,8,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    if (id.x >= Width || id.y >= Height)
        return;

    int2 p = int2(id.xy);
    int3 coord = int3(p, 0);

    // Obstacles: keep them fixed (no wave)
    float obs = Obstacles.Load(coord);
    if (obs > 0.5)
    {
        NState[p] = 0.0;
        return;
    }

    // Current and previous height
    float u     = State.Load(coord);
    float uPrev = PState.Load(coord);

    // Neighbour coords, clamped at borders
    int2 left   = int2(max(p.x - 1, 0),            p.y);
    int2 right  = int2(min(p.x + 1, (int)Width-1), p.y);
    int2 down   = int2(p.x, max(p.y - 1, 0));
    int2 up     = int2(p.x, min(p.y + 1, (int)Height-1));

    float uL = State.Load(int3(left, 0));
    float uR = State.Load(int3(right, 0));
    float uD = State.Load(int3(down, 0));
    float uU = State.Load(int3(up, 0));

    // Laplacian
    float lap = (uL + uR + uD + uU) - 4.0 * u;

    // Local wave speed from depth map
    float cLocal = GetLocalSpeed(coord);
    float c2dt2  = cLocal * cLocal * dt * dt;

    // Simple damped wave update:
    // u_next = (2 - D*dt) u - (1 - D*dt) uPrev + c^2 dt^2 lap
    float next = (2.0 - D * dt) * u - (1.0 - D * dt) * uPrev + c2dt2 * lap;

    NState[p] = next;
}

This is not a perfect discretization of the variable-c wave equation (that would use c at the cell faces in the divergence term), but it’s more than good enough visually and will give you the “bending toward shallow water” behavior.

3. Fixing / updating the Combine kernel

Let’s also fix Combine to be legal HLSL and compatible with the new Sources type:

[numthreads(8,8,1)]
void Combine (uint3 id : SV_DispatchThreadID)
{
    if (id.x >= Width || id.y >= Height)
        return;

    int3 coord = int3(id.xy, 0);

    float stateVal = State.Load(coord);
    float srcVal   = Sources.Load(coord); // now float instead of float4

    NState[id.xy] = stateVal + srcVal;
}

You can use this as a pass that:

Adds impulses from Sources into the state,
Or blends in any externally generated pattern.

4. Unity side: what you need to hook up

On the C# side (rough sketch):

int width = simWidth;
int height = simHeight;

waveCS.SetInt("Width", width);
waveCS.SetInt("Height", height);

waveCS.SetFloat("dt", dt);
waveCS.SetFloat("cShallow", cShallow);
waveCS.SetFloat("cDeep", cDeep);
waveCS.SetFloat("D", damping);

waveCS.SetTexture(csMainKernel, "PState", prevStateRT);
waveCS.SetTexture(csMainKernel, "State",  currStateRT);
waveCS.SetTexture(csMainKernel, "NState", nextStateRT);

waveCS.SetTexture(csMainKernel, "Obstacles", obstaclesRT);
waveCS.SetTexture(csMainKernel, "Depth",     depthRT);
waveCS.SetTexture(combineKernel, "Sources",  sourcesRT);

int groupsX = (width  + 7) / 8;
int groupsY = (height + 7) / 8;

waveCS.Dispatch(csMainKernel, groupsX, groupsY, 1);
// then Combine if you want source injection:
waveCS.Dispatch(combineKernel, groupsX, groupsY, 1);

Where:

depthRT is your depth map (aligned with the sim grid).
Values in depthRT are in [0,1]; 0 = shallow, 1 = deep (or flip it, but then adjust the lerp).

5. Why this gives you refraction

Because cLocal comes from Depth:

Deeper regions → larger cLocal → waves propagate faster.
Shallower regions → smaller cLocal → waves slow down.

When a wavefront hits the shallows at an angle, the part that reaches shallow depth first slows down first, while the rest keeps going fast for a bit, so the front rotates and lines up with the depth contours (and thus with the shore).

The cool part: you can literally paint your depth map and watch the waves “flow around” headlands, focus into bays, etc., just from the spatially varying c(x,y).

If you want to push this even further, next steps could be:

Use the gradient of the depth to add extra directional bias (simulate mild current / shoaling).
Make damping depend on depth too, so waves die more near the beach.
Add a simple nonlinearity near very shallow depths to mimic wave breaking (clamp or compress large amplitudes).

If you tell me how you’re generating / storing the depth map (e.g. from real geometry, hand-painted texture, etc.), I can tailor the mapping Depth → c(x,y) to give you stronger or subtler refraction.

Awesome, let’s turn this into a “feels-like-the-ocean” sim 😄
I’ll implement all 3 ideas directly in the compute shader and explain what each one is doing.

I’ll assume:

Depth is a Texture2D<float> in [0,1]:
0 = very shallow / beach, 1 = deep water.
You’re okay with a few more texture reads per step (still cheap enough for many use cases).

1. Full updated `CSMain` with all 3 effects

New parameters & helpers

#pragma kernel CSMain
#pragma kernel Combine

// States
Texture2D<float>   PState;     // previous state
Texture2D<float>   State;      // current state
RWTexture2D<float> NState;     // next state

// Inputs
Texture2D<float>   Sources;    // wave sources (scalar)
Texture2D<float>   Obstacles;  // 0 = water, 1 = solid
Texture2D<float>   Depth;      // 0 = shallow, 1 = deep

cbuffer WaveParams
{
    float dt;

    // Speed vs depth
    float cShallow;              // speed at Depth=0
    float cDeep;                 // speed at Depth=1

    // Damping vs depth
    float DShallow;              // damping at Depth=0 (near beach)
    float DDeep;                 // damping at Depth=1 (deep water)

    // Depth-gradient-based advection (shoaling / “current”)
    float DepthAdvectionStrength; // e.g. small value like 0.2

    // Breaking / nonlinearity near very shallow water
    float BreakDepth;            // depth threshold (0..1) below which breaking kicks in
    float BreakMaxAmplitude;     // approximate max wave height in shallow region
    float BreakSharpness;        // how aggressively to compress (e.g. 2-4)

    uint  Width;
    uint  Height;
};

float GetLocalSpeed(float depth01)
{
    // Depth in [0,1], map to speed
    return lerp(cShallow, cDeep, depth01);
}

float GetLocalDamping(float depth01)
{
    // More damping in shallow water (DShallow > DDeep usually)
    return lerp(DShallow, DDeep, depth01);
}

// Simple breaking model: compress large amplitudes near very shallow depths
float ApplyBreaking(float next, float depth01)
{
    // How shallow are we? 0 = deep (no breaking), 1 = at / above BreakDepth (max breaking)
    float shallowFactor = saturate((BreakDepth - depth01) / max(BreakDepth, 1e-4));

    if (shallowFactor <= 0.0)
        return next; // no breaking in deep water

    float amp = abs(next);
    if (amp <= BreakMaxAmplitude)
        return next; // under limit, do nothing

    // Soft-compression rather than hard clamp:
    // As amp grows above BreakMaxAmplitude, tanh squashes it.
    float t = amp / BreakMaxAmplitude;
    float compressedAmp = BreakMaxAmplitude * tanh(t * BreakSharpness);

    // Blend between original and compressed based on how shallow it is
    float finalAmp = lerp(amp, compressedAmp, shallowFactor);

    return (next >= 0.0 ? +finalAmp : -finalAmp);
}

Depth-aware, gradient-biased wave propagation

[numthreads(8,8,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
    if (id.x >= Width || id.y >= Height)
        return;

    int2 p = int2(id.xy);
    int3 coord = int3(p, 0);

    // Obstacles: keep them fixed (no wave)
    float obs = Obstacles.Load(coord);
    if (obs > 0.5)
    {
        NState[p] = 0.0;
        return;
    }

    // Current and previous height
    float u     = State.Load(coord);
    float uPrev = PState.Load(coord);

    // Neighbour coordinates (clamped at borders)
    int2 left   = int2(max(p.x - 1, 0),            p.y);
    int2 right  = int2(min(p.x + 1, (int)Width-1), p.y);
    int2 down   = int2(p.x, max(p.y - 1, 0));
    int2 up     = int2(p.x, min(p.y + 1, (int)Height-1));

    float uL = State.Load(int3(left, 0));
    float uR = State.Load(int3(right, 0));
    float uD = State.Load(int3(down, 0));
    float uU = State.Load(int3(up, 0));

    // Laplacian (4-neighbour)
    float lap = (uL + uR + uD + uU) - 4.0 * u;

    // Depth at center and neighbours
    float dC = Depth.Load(coord);
    float dL = Depth.Load(int3(left, 0));
    float dR = Depth.Load(int3(right, 0));
    float dD = Depth.Load(int3(down, 0));
    float dU = Depth.Load(int3(up, 0));

    // Local wave speed and damping from depth
    float cLocal = GetLocalSpeed(dC);
    float DLocal = GetLocalDamping(dC);

    float c2dt2 = cLocal * cLocal * dt * dt;

    // Basic damped wave update:
    // u_next = (2 - D*dt) u - (1 - D*dt) uPrev + c^2 dt^2 lap
    float damp1 = (2.0 - DLocal * dt);
    float damp2 = (1.0 - DLocal * dt);
    float next  = damp1 * u - damp2 * uPrev + c2dt2 * lap;

    // -------------------------------
    // Depth-gradient-based advection
    // -------------------------------

    // Gradient of depth (points toward deeper water if Depth increases)
    float2 gradDepth = float2((dR - dL) * 0.5, (dU - dD) * 0.5);

    // Gradient of height (wave slope)
    float2 gradU = float2((uR - uL) * 0.5, (uU - uD) * 0.5);

    // Make a "current" proportional to depth gradient
    // Positive DepthAdvectionStrength: energy tends to slide toward deeper water.
    float2 velocity = DepthAdvectionStrength * gradDepth;

    // Advection term: u_t + v·∇u = 0 => u_t = -v·∇u
    float advection = -dot(velocity, gradU);

    // Semi-hacky: just add a dt * advection contribution to the next step.
    next += dt * advection;

    // -------------------------------
    // Breaking / nonlinearity in shallow water
    // -------------------------------
    next = ApplyBreaking(next, dC);

    NState[p] = next;
}

Combine kernel (unchanged except for fixes)

[numthreads(8,8,1)]
void Combine (uint3 id : SV_DispatchThreadID)
{
    if (id.x >= Width || id.y >= Height)
        return;

    int3 coord = int3(id.xy, 0);

    float stateVal = State.Load(coord);
    float srcVal   = Sources.Load(coord);

    NState[id.xy] = stateVal + srcVal;
}

2. What each effect is doing (and how to tune it)

2.1. Depth gradient → directional bias (shoaling / “current”)

Core bits:

float2 gradDepth = float2((dR - dL) * 0.5, (dU - dD) * 0.5);
float2 gradU     = float2((uR - uL) * 0.5, (uU - uD) * 0.5);
float2 velocity  = DepthAdvectionStrength * gradDepth;
float  advection = -dot(velocity, gradU);
next += dt * advection;

gradDepth tells you which way the depth is increasing.
- If Depth is 0 in shallow and 1 in deep:
  - gradDepth points from shallow → deep.
We pretend there’s a “current” velocity v = DepthAdvectionStrength * gradDepth.
- So flow tends to follow the slope toward deeper water.
gradU is the slope of the wave surface.
For a simple advection equation:
$ u_t + \mathbf{v} \cdot \nabla u = 0 \Rightarrow u_t = -\mathbf{v} \cdot \nabla u $.
We add dt * advection to approximate that.

Effect visually:

Waves slightly lean / drift in the direction of deeper water.
In combination with variable-speed refraction, you get more “shoaling” behavior:
- crests bend and slide along contours instead of just straight-in refraction.

Tuning:

Start with DepthAdvectionStrength ≈ 0.1–0.3.
- 0 means no extra bias (only refraction from variable c).
- Larger values → more noticeable flow, can become “smearing” if too high.

2.2. Depth-dependent damping (waves die out near the beach)

Core bits:

float GetLocalDamping(float depth01)
{
    return lerp(DShallow, DDeep, depth01);
}

// ...

float DLocal = GetLocalDamping(dC);

float damp1 = (2.0 - DLocal * dt);
float damp2 = (1.0 - DLocal * dt);
float next  = damp1 * u - damp2 * uPrev + c2dt2 * lap;

DShallow = damping at Depth=0 (beach).
DDeep = damping at Depth=1 (deep water).
If you choose DShallow > DDeep, then:
- Waves lose more energy when they get into very shallow water.
- Out at sea, they’re almost undamped (or just a tiny damping for numerical safety).

Effect visually:

Waves “fade out” as they hit the shore instead of endlessly bouncing around.
Combined with breaking (next section), you get that classic dying, foamy shore zone feeling (minus the actual foam rendering).

Tuning:

Example starting values:
- DDeep ≈ 0.2
- DShallow ≈ 1.0 or higher
If the whole sim feels too lossy:
- Reduce both, but keep DShallow > DDeep.

2.3. Breaking / nonlinearity in very shallow water

Core bits:

float ApplyBreaking(float next, float depth01)
{
    float shallowFactor = saturate((BreakDepth - depth01) / max(BreakDepth, 1e-4));

    if (shallowFactor <= 0.0)
        return next; // deep water, no breaking

    float amp = abs(next);
    if (amp <= BreakMaxAmplitude)
        return next; // small waves, no breaking

    float t = amp / BreakMaxAmplitude;
    float compressedAmp = BreakMaxAmplitude * tanh(t * BreakSharpness);

    float finalAmp = lerp(amp, compressedAmp, shallowFactor);

    return (next >= 0.0 ? +finalAmp : -finalAmp);
}

Breaking only happens when:
- We’re shallower than BreakDepth, and
- Amplitude exceeds BreakMaxAmplitude.
shallowFactor:
- 0 in deep water,
- 1 at depth01 = 0 (very shallow).
- Allows smooth fade-in of breaking effect as you get close to shore.
tanh gives a soft clamp:
- For small t, tanh(t) ≈ t → no change.
- For big t, tanh(t) approaches 1 → amplitude asymptotes around BreakMaxAmplitude.
Then we blend between original amp and compressed amp based on shallowFactor:
- Deep water: no compression at all.
- Very shallow: fully compressed when above the threshold.

Effect visually:

In deep water: waves behave linearly and can get larger.
As they enter very shallow regions:
- Their height growth is limited, so tall crests start to “flatten” instead of just growing.
- This feels like crude wave breaking without a full fluid sim.

Tuning:

BreakDepth: where in the depth range you want breaking to start.
- If Depth is 0..1 with 0 = beach:
  - BreakDepth ≈ 0.2–0.4 means: start breaking in last 20–40% of shallows.
BreakMaxAmplitude: max height you’re OK with near shore.
- Adjust until it “feels right” vs your height scaling.
BreakSharpness:
- 1 → very gentle compression.
- 3-4 → more “hard” breaking/clipping look.

If you want even more cartoon-ish breaking, you can replace the tanh compression with a straight clamp:

float compressedAmp = min(amp, BreakMaxAmplitude);

But tanh keeps it smoother and more natural.

3. What to do on the C# side

You’ll need to set the new uniforms:

waveCS.SetFloat("dt", dt);

// speed vs depth
waveCS.SetFloat("cShallow", cShallow);
waveCS.SetFloat("cDeep",    cDeep);

// damping vs depth
waveCS.SetFloat("DShallow", DShallow);
waveCS.SetFloat("DDeep",    DDeep);

// advection strength
waveCS.SetFloat("DepthAdvectionStrength", depthAdvectionStrength);

// breaking
waveCS.SetFloat("BreakDepth",        breakDepth);
waveCS.SetFloat("BreakMaxAmplitude", breakMaxAmplitude);
waveCS.SetFloat("BreakSharpness",    breakSharpness);

// size
waveCS.SetInt("Width",  width);
waveCS.SetInt("Height", height);

Then bind your textures like before (State, PState, NState, Depth, Obstacles, Sources) and dispatch.

If you want, next we can go one level nerdier and tweak the update so the variable c(x,y) is treated in a more “correct PDE” way (c at cell faces, etc.) for better stability and less anisotropy—but honestly, for visual waves this combo (refraction from depth, gradient-driven drift, depth-dependent damping, breaking) is already a big, juicy upgrade.