環境

Unity2022.2.14f1

概要

炎のシミュレーションのテストです。

参考の動画を見つつ、ナビエストークスを「poisson separable filter」を使用して行ってみました。正直正しいのかわからないです。

ヤコビ反復32回分の代替になると解説している気がする。

コード

using UnityEngine;
using UnityEngine.UI;

public class TestFluid : MonoBehaviour
{
    public class DoubleRenderTexture
    {
        private RenderTexture[] _textures = new RenderTexture[2];
        public RenderTexture current => _textures[0];
        public RenderTexture back => _textures[1];
        private int _width = 0;
        private int _height = 0;
        public int width => _width;
        public int height => _height;
        public DoubleRenderTexture(int width, int height)
        {
            _width = width;
            _height = height;
            for (int i = 0; i < _textures.Length; i++)
            {
                _textures[i] = new RenderTexture(width, height, 0, RenderTextureFormat.ARGBHalf, RenderTextureReadWrite.Linear);
                _textures[i].enableRandomWrite = true;
                _textures[i].Create();
            }
        }

        public void Flip()
        {
            RenderTexture temp = _textures[0];
            _textures[0] = _textures[1];
            _textures[1] = temp;
        }

        public void Release()
        {
            foreach (var texture in _textures)
                texture.Release();
        }
    }

    [SerializeField] Material _material = null;
    [SerializeField] private ComputeShader _shader = null;
    [SerializeField] private RawImage _densityPreview = null;
    [SerializeField] private RawImage _velocityPreview = null;
    [SerializeField] private float _attenuation = 0.99f;
    [SerializeField] private float _deisityAttenuation = 0.985f;
    [SerializeField] private Vector2Int _vectorNum = new Vector2Int(16, 16);
    [SerializeField] private Vector2Int _densityNum = new Vector2Int(128, 128);

    private RenderTexture _divergenceTexture = null;
    private DoubleRenderTexture _velocityBuffer = null;
    private DoubleRenderTexture _pressureBuffer = null;
    private DoubleRenderTexture _densityBuffer = null;

    private int _knUpdateAdvection;
    private int _knInteractionForce;
    private int _knInteractionDensityForce;
    private int _knUpdateDivergence;
    private int _knUpdatePressurePsfH;
    private int _knUpdatePressurePsfV;
    private int _knUpdateVelocity;
    private int _knUpdateDensity;
    private static int _spTime = Shader.PropertyToID("_Time");
    private static int _spDeltaTime = Shader.PropertyToID("_DeltaTime");
    private static int _spSourceVelocity = Shader.PropertyToID("_SourceVelocity");
    private static int _spResultVelocity = Shader.PropertyToID("_ResultVelocity");
    private static int _spSourcePressure = Shader.PropertyToID("_SourcePressure");
    private static int _spResultPressure = Shader.PropertyToID("_ResultPressure");
    private static int _spResultDivergence = Shader.PropertyToID("_ResultDivergence");
    private static int _spSourceDensity = Shader.PropertyToID("_SourceDensity");
    private static int _spResultDensity = Shader.PropertyToID("_ResultDensity");
    private static int _spSimWidth = Shader.PropertyToID("_SimWidth");
    private static int _spSimHeight = Shader.PropertyToID("_SimHeight");
    private static int _spSimForceStart = Shader.PropertyToID("_SimForceStart");
    private static int _spDensityForceStart = Shader.PropertyToID("_DensityForceStart");
    private static int _spAttenuation = Shader.PropertyToID("_Attenuation");
    private static int _spDeisityAttenuation = Shader.PropertyToID("_DeisityAttenuation");
    private static int _spDensityWidth = Shader.PropertyToID("_DensityWidth");
    private static int _spDensityHeight = Shader.PropertyToID("_DensityHeight");
    private static int _spVectorFieldTex = Shader.PropertyToID("_VectorFieldTex");
    private static int _spDensityTex = Shader.PropertyToID("_DensityTex");

    private void Start()
    {
        Setup();
    }

    private void OnValidate()
    {
        _vectorNum.x = Mathf.Max(_vectorNum.x, 8);
        _vectorNum.y = Mathf.Max(_vectorNum.y, 8);
        _densityNum.x = Mathf.Max(_densityNum.x, 8);
        _densityNum.y = Mathf.Max(_densityNum.y, 8);

        DestroyBuffers();
        Setup();
    }

    private void OnDestroy()
    {
        DestroyBuffers();
    }

    private void Setup()
    {
        CreateBuffers();

        _knUpdateAdvection = _shader.FindKernel("UpdateAdvection");
        _knInteractionForce = _shader.FindKernel("InteractionForce");
        _knInteractionDensityForce = _shader.FindKernel("InteractionDensityForce");
        _knUpdateDivergence = _shader.FindKernel("UpdateDivergence");
        _knUpdatePressurePsfH = _shader.FindKernel("UpdatePressurePsfH");
        _knUpdatePressurePsfV = _shader.FindKernel("UpdatePressurePsfV");
        _knUpdateVelocity = _shader.FindKernel("UpdateVelocity");
        _knUpdateDensity = _shader.FindKernel("UpdateDensity");
    }

    private void CreateBuffers()
    {
        _velocityBuffer = new DoubleRenderTexture(_vectorNum.x, _vectorNum.y);
        _pressureBuffer = new DoubleRenderTexture(_vectorNum.x, _vectorNum.y);
        _densityBuffer = new DoubleRenderTexture(_densityNum.x, _densityNum.y);
        _divergenceTexture = new RenderTexture(_vectorNum.x, _vectorNum.y, 0, RenderTextureFormat.ARGBHalf, RenderTextureReadWrite.Linear);
        _divergenceTexture.enableRandomWrite = true;
        _divergenceTexture.Create();
    }

    private void DestroyBuffers()
    {
        if(_divergenceTexture != null)
            _divergenceTexture.Release();
        if(_velocityBuffer != null)
            _velocityBuffer.Release();
        if(_pressureBuffer != null)
            _pressureBuffer.Release();
        if(_densityBuffer != null)
            _densityBuffer.Release();
    }

    private void Update()
    {
        _shader.SetFloat(_spSimWidth, _vectorNum.x);
        _shader.SetFloat(_spSimHeight, _vectorNum.y);
        _shader.SetFloat(_spDeltaTime, Time.deltaTime);
        _shader.SetFloat(_spAttenuation, _attenuation);
        _shader.SetFloat(_spDeisityAttenuation, _deisityAttenuation);
        _shader.SetFloat(_spTime, Time.time);
        _shader.SetFloat(_spDensityWidth, _densityNum.x);
        _shader.SetFloat(_spDensityHeight, _densityNum.y);

        UpdateAdvection();
        InteractionForce();
        UpdateDivergence();
        UpdatePressure();
        UpdateVelocity();
        UpdateDensity();

        _velocityPreview.texture = _velocityBuffer.current;
        _densityPreview.texture = _densityBuffer.current;
        _material.SetTexture(_spVectorFieldTex, _velocityBuffer.current);
        _material.SetTexture(_spDensityTex, _densityBuffer.current);
    }

    // 移流
    private void UpdateAdvection()
    {
        _shader.SetTexture(_knUpdateAdvection, _spSourceVelocity, _velocityBuffer.current);
        _shader.SetTexture(_knUpdateAdvection, _spResultVelocity, _velocityBuffer.back);
        _shader.Dispatch(_knUpdateAdvection, _velocityBuffer.width / 8, _velocityBuffer.height / 8, 1);
        _velocityBuffer.Flip();
    }

    // 外力
    private void InteractionForce()
    {
        var width = _vectorNum.x / 8;
        var halfWidth = width / 2;
        var halfX = _vectorNum.x / 2;
        var start = halfX - halfWidth;
        _shader.SetInt(_spSimForceStart, start);
        _shader.SetTexture(_knInteractionForce, _spResultVelocity, _velocityBuffer.current);
        _shader.Dispatch(_knInteractionForce, width, 1, 1);

        // Density
        width = _densityNum.x / 8;
        halfWidth = width / 2;
        halfX = _densityNum.x / 2;
        start = halfX - halfWidth;
        _shader.SetInt(_spDensityForceStart, start);
        _shader.SetTexture(_knInteractionDensityForce, _spResultDensity, _densityBuffer.current);
        _shader.Dispatch(_knInteractionDensityForce, width, 1, 1);
    }

    // 発散
    private void UpdateDivergence()
    {
        _shader.SetTexture(_knUpdateDivergence, _spSourceVelocity, _velocityBuffer.current);
        _shader.SetTexture(_knUpdateDivergence, _spResultDivergence, _divergenceTexture);
        _shader.Dispatch(_knUpdateDivergence, _divergenceTexture.width / 8, _divergenceTexture.height / 8, 1);
    }

    // 圧力
    private void UpdatePressure()
    {
        // Horizontal
        _shader.SetTexture(_knUpdatePressurePsfH, _spSourcePressure, _pressureBuffer.current);
        _shader.SetTexture(_knUpdatePressurePsfH, _spResultPressure, _pressureBuffer.back);
        _shader.SetTexture(_knUpdatePressurePsfH, _spResultDivergence, _divergenceTexture);
        _shader.Dispatch(_knUpdatePressurePsfH, _pressureBuffer.width / 8, _pressureBuffer.height / 8, 1);
        // Vertical
        _shader.SetTexture(_knUpdatePressurePsfV, _spSourcePressure, _pressureBuffer.current);
        _shader.SetTexture(_knUpdatePressurePsfV, _spResultPressure, _pressureBuffer.back);
        _shader.SetTexture(_knUpdatePressurePsfV, _spResultDivergence, _divergenceTexture);
        _shader.Dispatch(_knUpdatePressurePsfV, _pressureBuffer.width / 8, _pressureBuffer.height / 8, 1);

        _pressureBuffer.Flip();
    }

    // 速度
    private void UpdateVelocity()
    {
        _shader.SetTexture(_knUpdateVelocity, _spSourceVelocity, _velocityBuffer.current);
        _shader.SetTexture(_knUpdateVelocity, _spSourcePressure, _pressureBuffer.current);
        _shader.SetTexture(_knUpdateVelocity, _spResultVelocity, _velocityBuffer.back);
        _shader.Dispatch(_knUpdateVelocity, _velocityBuffer.width / 8, _velocityBuffer.height / 8, 1);
        _velocityBuffer.Flip();
    }

    // 密度
    private void UpdateDensity()
    {
        _shader.SetTexture(_knUpdateDensity, _spSourceDensity, _densityBuffer.current);
        _shader.SetTexture(_knUpdateDensity, _spSourceVelocity, _velocityBuffer.current);
        _shader.SetTexture(_knUpdateDensity, _spResultDensity, _densityBuffer.back);
        _shader.Dispatch(_knUpdateDensity, _densityBuffer.width / 8, _densityBuffer.height / 8, 1);
        _densityBuffer.Flip();
    }
}

#pragma kernel UpdateAdvection
#pragma kernel InteractionForce
#pragma kernel InteractionDensityForce
#pragma kernel UpdateDivergence
#pragma kernel UpdatePressurePsfH
#pragma kernel UpdatePressurePsfV
#pragma kernel UpdateVelocity
#pragma kernel UpdateDensity

Texture2D<float4> _SourceVelocity;
Texture2D<float4> _SourcePressure;
Texture2D<float4> _SourceDensity;

RWTexture2D<float4> _ResultVelocity;
RWTexture2D<float4> _ResultPressure;
RWTexture2D<float4> _ResultDivergence;
RWTexture2D<float4> _ResultDensity;

#define SAMPLE _LinearClamp
SamplerState _LinearClamp;
SamplerState _PointClamp;

float _Time;
float _DeltaTime;
float _SimWidth;
float _SimHeight;
uint _SimForceStart;
uint _DensityForceStart;
float _DensityWidth;
float _DensityHeight;
float _Attenuation;
float _DeisityAttenuation;

#define PI 3.1415926538
#define PI2 (PI * 2)

float RandomRange(float2 Seed, float Min, float Max)
{
    float randomno =  frac(sin(dot(Seed, float2(12.9898, 78.233)))*43758.5453);
    return lerp(Min, Max, randomno);
}

// 移流
[numthreads(8,8,1)]
void UpdateAdvection(uint2 id : SV_DispatchThreadID)
{
    float w = _SimWidth;
    float h = _SimHeight;

    float3 px = float3(1.0/w, 1.0/h, 0.0);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

    float2 velocity = _SourceVelocity.SampleLevel(SAMPLE, uv, 0).xy;
    float2 result = _SourceVelocity.SampleLevel(SAMPLE, uv - velocity * _DeltaTime, 0).xy;

    _ResultVelocity[id] = float4(result, 0.0, 1.0);
}

// 外力
[numthreads(1,1,1)]
void InteractionForce(uint2 id : SV_DispatchThreadID)
{
    id.x += _SimForceStart;
    const float power = 1.5;
    float3 vec = float3(0.0, power, 0.0);
    _ResultVelocity[id] = float4(vec, 1.0);

    // 適当な風力
    const float speed = 5.0;
    float rad = _Time * speed;
    id.x -= _SimForceStart;
    id.y = _SimHeight / 2;
    vec.xy = float2(cos(rad), sin(rad)) * power * 1.0;
    _ResultVelocity[id] = float4(vec, 1.0);
}

[numthreads(1,1,1)]
void InteractionDensityForce(uint2 id : SV_DispatchThreadID)
{
    id.x += _DensityForceStart;
    _ResultDensity[id] = RandomRange(_Time + id.x, 0.3, 1.0);
}

// 発散
[numthreads(8,8,1)]
void UpdateDivergence(uint2 id : SV_DispatchThreadID)
{
    float w = _SimWidth;
    float h = _SimHeight;

    float3 px = float3(1.0 / w, 1.0 / h, 0);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

    float x0 = _SourceVelocity.SampleLevel(SAMPLE, uv - px.xz, 0).x;
    float x1 = _SourceVelocity.SampleLevel(SAMPLE, uv + px.xz, 0).x;
    float y0 = _SourceVelocity.SampleLevel(SAMPLE, uv - px.zy, 0).y;
    float y1 = _SourceVelocity.SampleLevel(SAMPLE, uv + px.zy, 0).y;

    float divergence = (x1 - x0 + y1 - y0);

    _ResultDivergence[id] = float4(divergence.xx, 0.0, 1.0);
}

static const float poissonFilter[7] = {
    .57843719174,
    .36519596949,
    .23187988879,
    .14529589353,
    .08816487385,
    .05184872885,
    .02906462467
};

// 圧力 H
[numthreads(8,8,1)]
void UpdatePressurePsfH(uint2 id : SV_DispatchThreadID)
{
    float w = _SimWidth;
    float h = _SimHeight;

    float3 px = float3(1.0 / w, 1.0 / h, 0);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;
    int i;

    float p0 = 0.0;
    [unroll]
    for (i = -6; i <= 6; i++) {
        float x = _SourcePressure.SampleLevel(_LinearClamp, uv + float2(px.x * i, 0), 0).r;
        p0 += poissonFilter[abs(i)] * x;
    }
    _ResultPressure[id] = float4(p0, 0.0, 0.0, 0.0);
}

// 圧力 V
[numthreads(8,8,1)]
void UpdatePressurePsfV(uint2 id : SV_DispatchThreadID)
{
    float w = _SimWidth;
    float h = _SimHeight;

    float3 px = float3(1.0 / w, 1.0 / h, 0);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;
    int i;

    float p1 = 0.0;
    [unroll]
    for (i = -6; i <= 6; i++) {
        float y = _SourcePressure.SampleLevel(_LinearClamp, uv + float2(0, px.y * i), 0).r;
        p1 += poissonFilter[abs(i)] * y;
    }

    float p0 = _ResultPressure[id].x;
    float d = _ResultDivergence[id].r;
    const float scale = 0.5;
    float relaxed = (p1 - d) * 0.25 * scale;

    _ResultPressure[id] = float4(relaxed.xx, 0.0, 1.0);
}

// 速度
[numthreads(8,8,1)]
void UpdateVelocity(uint2 id : SV_DispatchThreadID)
{
    float w = _SimWidth;
    float h = _SimHeight;

    float3 px = float3(1.0 / w, 1.0 / h, 0);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

    float x0 = _SourcePressure.SampleLevel(SAMPLE, uv - px.xz, 0).r;
    float x1 = _SourcePressure.SampleLevel(SAMPLE, uv + px.xz, 0).r;
    float y0 = _SourcePressure.SampleLevel(SAMPLE, uv - px.zy, 0).r;
    float y1 = _SourcePressure.SampleLevel(SAMPLE, uv + px.zy, 0).r;

    float2 v = _SourceVelocity.SampleLevel(SAMPLE, uv, 0).xy;
    float4 v2 = float4((v - (float2(x1, y1) - float2(x0, y0)) * 0.5), 0.0, 1.0);
    v2 *= _Attenuation;

    _ResultVelocity[id] = v2;
}

// 密度
[numthreads(8,8,1)]
void UpdateDensity(uint2 id : SV_DispatchThreadID)
{
    float dw = _DensityWidth;
    float dh = _DensityHeight;
    float2 baseUv = float2(id.x / dw, id.y / dh);
    float3 dpx = float3(1.0 / dw, 1.0 / dh, 0);
    float2 duv = baseUv + dpx.xy * 0.5;

    float vw = _SimWidth;
    float vh = _SimHeight;
    float3 vpx = float3(1.0 / vw, 1.0 / vh, 0);
    const float scale = 0.075;
    float2 vuv = baseUv + vpx.xy * scale;

    float2 vel = _SourceVelocity.SampleLevel(SAMPLE, vuv, 0).xy;

    float4 col = _SourceDensity.SampleLevel(SAMPLE, duv - vel * _DeltaTime, 0);
    col *= _DeisityAttenuation;
    _ResultDensity[id] = col;
}

Shader "Custom/Fire"
{
    Properties
    {
        _GradientTex ("Gradient Texture", 2D) = "white" {}
    }

    CGINCLUDE
    #include "UnityCG.cginc"

    struct appdata
    {
        float4 vertex : POSITION;
        float2 uv : TEXCOORD0;
    };

    struct v2f
    {
        float4 vertex : SV_POSITION;
        float2 uv : TEXCOORD0;
    };

    sampler2D _VectorFieldTex;
    float4 _VectorFieldTex_TexelSize;
    sampler2D _DensityTex;
    float4 _DensityTex_TexelSize;
    sampler2D _GradientTex;
 
    static const float PI = 3.1415927;
    static const int   ARROW_V_STYLE = 1;
    static const int   ARROW_LINE_STYLE = 2;
    static const int   ARROW_STYLE = ARROW_LINE_STYLE;
    static const float ARROW_HEAD_ANGLE = 45.0 * PI / 180.0;
    static const float ARROW_SHAFT_THICKNESS = 3.0;

    static const float ARROW_BASE_SIZE = 1024.0;
    #define ARROW_TILE_SIZE (ARROW_BASE_SIZE / (float)_VectorFieldTex_TexelSize.z)
    #define ARROW_HEAD_LENGTH (ARROW_TILE_SIZE / 6.0)

 
    float2 ArrowTileCenterCoord(float2 pos)
    {
        return (floor(pos / ARROW_TILE_SIZE) + 0.5) * ARROW_TILE_SIZE;
    }
 
    float Arrow(float2 p, float2 v)
    {
        p -= ArrowTileCenterCoord(p);
        float mag_v = length(v), mag_p = length(p);  
        if (mag_v > 0.0)
        {
            float2 dir_p = p / mag_p, dir_v = v / mag_v;
            mag_v = clamp(mag_v, 5.0, ARROW_TILE_SIZE / 2.0);
            v = dir_v * mag_v;
            float dist;
            if (ARROW_STYLE == ARROW_LINE_STYLE)
            {
                dist = max(ARROW_SHAFT_THICKNESS / 4.0 - max(abs(dot(p, float2(dir_v.y, -dir_v.x))),
                    abs(dot(p, dir_v)) - mag_v + ARROW_HEAD_LENGTH / 2.0),
                    min(0.0, dot(v - p, dir_v) - cos(ARROW_HEAD_ANGLE / 2.0) * length(v - p)) * 2.0 +
                    min(0.0, dot(p, dir_v) + ARROW_HEAD_LENGTH - mag_v));
            }
            else
            {
                dist = min(0.0, mag_v - mag_p) * 2.0 +
                        min(0.0, dot(normalize(v - p), dir_v) - cos(ARROW_HEAD_ANGLE / 2.0)) * 2.0 * length(v - p) +
                        min(0.0, dot(p, dir_v) + 1.0) +
                        min(0.0, cos(ARROW_HEAD_ANGLE / 2.0) - dot(normalize(v * 0.33 - p), dir_v)) * mag_v * 0.8;
            }
            return clamp(1.0 + dist, 0.0, 1.0);
        }
        else
        {
            return max(0.0, 1.2 - mag_p);
        }
    }
 
    float2 field(float2 pos)
    {
        return tex2D(_VectorFieldTex, pos / ARROW_BASE_SIZE);
    }

    ENDCG

    SubShader
    {
        Tags
        { 
            "RenderType" = "Transparent"
            "Queue" = "Transparent"
        }

        ZWrite Off
        Lighting Off
        Blend SrcAlpha OneMinusSrcAlpha
        Pass
        {
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            v2f vert(appdata v)
            {
                v2f o;
                o.vertex = UnityObjectToClipPos(v.vertex);
                o.uv = v.uv;
                return o;
            }
 
            float4 frag(v2f i) : SV_Target
            {
                float2 fragCoord = i.uv * float2(ARROW_BASE_SIZE,ARROW_BASE_SIZE);

                float4 col;

                float arrow = (1.0 - Arrow(fragCoord.xy, field(ArrowTileCenterCoord(fragCoord.xy)) * ARROW_TILE_SIZE * 0.4));
                float4 densityCol = tex2D(_GradientTex, float2(saturate(tex2D(_DensityTex, i.uv).r), 0.5));
                col.rgb = densityCol.rgb * arrow;
                col.a = max(densityCol.a, 1.0 - arrow);
                return col;
            }
            ENDCG
        }
    }
}