LYGIA Shader Library

Dependencies:

• lygia/v1.1.6/math/saturate.glsl
• lygia/v1.1.6/sample.glsl

Use:

<vec2> fluidSolver(<SAMPLER_TYPE> tex, <vec2> st, <vec2> pixel, <vec2> force)

Check it on Github

#ifndef FLUIDSOLVER_DT
#define FLUIDSOLVER_DT 0.15
#endif

#ifndef FLUIDSOLVER_DX
#define FLUIDSOLVER_DX 1.0
#endif

// higher this threshold, lower the viscosity (max .8)
#ifndef FLUIDSOLVER_VISCOSITY
#define FLUIDSOLVER_VISCOSITY .16
#endif

// #ifndef FLUIDSOLVER_VELOCITY_DECAY
// #define FLUIDSOLVER_VELOCITY_DECAY 5e-6
// #endif 

#ifndef FLUIDSOLVER_SAMPLER_FNC
#define FLUIDSOLVER_SAMPLER_FNC(TEX, UV) SAMPLER_FNC(TEX, UV)
#endif

#ifndef FNC_FLUIDSOLVER
#define FNC_FLUIDSOLVER

vec4 fluidSolver(SAMPLER_TYPE tex, vec2 st, vec2 pixel, vec2 force) {
    const float k = .2;
    const float s = k/FLUIDSOLVER_DT;
    const float dx = FLUIDSOLVER_DX;

    // Data
    vec4 d = FLUIDSOLVER_SAMPLER_FNC(tex, st);
    vec4 dR = FLUIDSOLVER_SAMPLER_FNC(tex, st + vec2(pixel.x, 0.));
    vec4 dL = FLUIDSOLVER_SAMPLER_FNC(tex, st - vec2(pixel.x, 0.));
    vec4 dT = FLUIDSOLVER_SAMPLER_FNC(tex, st + vec2(0., pixel.y));
    vec4 dB = FLUIDSOLVER_SAMPLER_FNC(tex, st - vec2(0., pixel.y));

    // Delta Data
    vec4 ddx = (dR - dL).xyzw; // delta data on X
    vec4 ddy = (dT - dB).xyzw; // delta data on Y
    float divergence = (ddx.x + ddy.y) * 0.5;
    divergence = (ddx.x + ddy.y) / (2.0 * dx * dx);

    // Solving for density with one jacobi iteration 
    float a = 1.0 / (dx * dx);
    d.z = 1.0 / ( -4.0 * a ) * ( divergence - a * (dT.z + dR.z + dB.z + dL.z));

    #ifdef FLUIDSOLVER_VISCOSITY
    // Solving for velocity
    vec2 laplacian = dR.xy + dL.xy + dT.xy + dB.xy - 4.0 * d.xy;
    vec2 viscosityForce = FLUIDSOLVER_VISCOSITY * laplacian;
    // Semi-lagrangian advection
    vec2 densityInvariance = s * vec2(ddx.z, ddy.z);
    vec2 was = st - FLUIDSOLVER_DT * d.xy * pixel;
    d.xyw = FLUIDSOLVER_SAMPLER_FNC(tex, was).xyw;
    d.xy += FLUIDSOLVER_DT * (viscosityForce - densityInvariance + force);
    #endif

    #ifdef FLUIDSOLVER_VELOCITY_DECAY
    d.xy = max(vec2(0.), abs(d.xy) - FLUIDSOLVER_VELOCITY_DECAY) * sign(d.xy);
    #endif

    // Vorticity confinement
    #ifdef FLUIDSOLVER_VORTICITY
    d.w = (dB.x - dT.x + dR.y - dL.y); // curl stored in the w channel
    vec2 vorticity = vec2(abs(dT.w) - abs(dB.w), abs(dL.w) - abs(dR.w));
    vorticity *= FLUIDSOLVER_VORTICITY/(length(vorticity) + 1e-5) * d.w;
    d.xy += vorticity;
    #endif

    // Boundary conditions
    d.xy *= smoothstep(0.5, 0.49,abs(st - 0.5));

    // Pack XY, Z and W data
    d.xy = clamp(d.xy, -0.999, 0.999);// * 0.5 + 0.5;
    d.zw = saturate(d.zw);
    return d;
}

#endif