Vala プログラミング

WebGPU プログラミング

おなが@京都先端科学大

Julia MeshCat SPH 流体シミュレーション (PointCloud)

Julia

Julia MeshCat ライブラリを用いて、水柱崩壊のシミュレーションを行ってみました。
流体のシミュレーションには、SPH (Smoothed Particle Hydrodynamics) 法を
使用しています。
f:id:onagat12:20200604174440g:plain

  • MeshCat ( Julia bindings to the MeshCat WebGL viewer)

　https://github.com/rdeits/MeshCat.jl
　MeshCat は Three.js を用いた3D Wiewer です。簡単にmeshの描画ができるので、
　MeshCat を利用しました。

  • SPH

　SPHによる流体シミュレーションは、様々なサイトで紹介されています。
　ここで使用した数式と数値は、以下の２つのサイトを参考にしています。
　1. SPH法による流体シミュレーション
　　https://soysoftware.sakura.ne.jp/archives/1559
　2. 粒子法のプログラム
　　http://kamonama.blogspot.com/2009/02/blog-post_23.html

  • 粒子の描画

　粒子はPointCloudで描画しています。
　PointCluodの描画

color = RGBA{Float32}(0.0, 1.0, 0.0, 1.0)
material = PointsMaterial(color = color, size = 0.05, vertexColors = 0)
setobject!(vis, PointCloud(positions), material)

　PointCloudの色の設定、ポイントサイズの設定は、PointMaterialで行います。

  • 画像のキャプチャー

　MeshCat ウィンドウのControls から、save_image を使ってキャプチャーします。
　Base.prompt()文（キー入力待ち）を挿入し、処理を停止させて、 １コマ１コマ
　キャプチャーしています。動画は200コマを使用。

プログラム

# MeshCat-sph.jl

using MeshCat
vis = Visualizer()
#open(vis)

using Blink
AtomShell.isinstalled() || AtomShell.install()
open(vis, Window())

using ColorTypes
using CoordinateTransformations
using LinearAlgebra
using GeometryTypes

mutable struct Particle
    position::Point3
    velocity::Point3
    density
    pressure
    force::Point3
end

# weight function
w_rho(r, h) = r < h ? (h^2 - r^2)^3 : 0.0
w_pressure(r,h) = r < h ? (h - r)^2 : 0.0
w_viscosity(r,h) = r < h ? (h- r) : 0.0

# Gravity
accel_g = [0.0, 0.0, -9.8]

# constants
mass = 0.00020543
pressure_stiffness = 1.0
rest_density = 600.0
h = 0.01
l0 = (mass / rest_density)^(1.0 / 3.0) # l0 = 0.0069958
simscale = 0.004
d = l0 / simscale * 0.86 # d = 1.5216
visc = 0.2
limit = 200.0
radius = 0.004
epsiron = 0.00001
extstiff = 10000.0
extdamp = 256.0

w_rho_coef = 315.0 / (64.0 * pi * h^9)
w_pressure_coef = 45.0 / (pi * h^6)
w_visc_coef = w_pressure_coef

init_min = [ 0.0,  0.0, 0.0]
pos_min = [0.0, 0.0, 0.0]
pos_max = [10.0, 25.0, 20.0] 

nx = 5; ny = 5; nz = 10

particle_list = []

for iz = 0:nz+1, iy = 0:ny+1, ix = 0:nx+1
    x = init_min[1] + d*ix; y = init_min[2] + d*iy; z = init_min[3] + d*iz

    position = Point3(x, y, z)
    velocity = Point3(0, 0, 0)
    density = 0.0
    pressure = 0.0
    force = Point3(0.0, 0.0, 0.0)
    push!(particle_list,
       Particle(position, velocity, density, pressure, force))
end
particles = length(particle_list)
   println("number of particles: ", particles)

function Simulation(DT)
   # neghiber_list
    neighbor_list = []
    for i = 1:particles
        position_i = particle_list[i].position
        neighbors = []
        for j = 1:particles
            j == i && continue
            position_j = particle_list[j].position
            distance = norm(position_j - position_i) * simscale
            if distance < h
                push!(neighbors, j)
            end
        end
        push!(neighbor_list, neighbors)
        #println("neighbors: ", neighbors)
    end

    # density(rho)
    for i = 1:particles
        position_i = particle_list[i].position

        sum = 0.0
        map(neighbor_list[i]) do j
            position_j = particle_list[j].position
            distance = norm(position_j - position_i) * simscale
            sum += w_rho(distance, h)
        end
        # density
        particle_list[i].density = sum * mass * w_rho_coef
        #println("density: ", i, " ", particle_list[i].density)

        # pressure
        particle_list[i].pressure =
            pressure_stiffness * (particle_list[i].density - rest_density)
        #println("pressure: ", particle_list[i].pressure)
    end
    
    # force: pressure gradient, viscosity and adj(force at wall)
    for i = 1:particles
        pressure_gradient = [0.0, 0.0, 0.0]
        visc_term = [0.0, 0.0, 0.0]
        position_i = particle_list[i].position
        
        map(neighbor_list[i]) do j
            position_j = particle_list[j].position
            distance = norm(position_j - position_i) * simscale
            #
            diffji =  (position_j - position_i) * simscale
            pweight = w_pressure_coef * w_pressure(distance, h) * diffji / distance
            
            pji = particle_list[j].pressure + particle_list[i].pressure

            rho_i = particle_list[i].density
            rho_j = particle_list[j].density
            pgradient = -0.5 * pji / (rho_i * rho_j) * pweight
            pressure_gradient .+= pgradient

            vweight = w_visc_coef * w_viscosity(distance, h)
            vji = particle_list[j].velocity - particle_list[i].velocity
            vterm = visc  * vji / (rho_i * rho_j) * vweight 
            
            visc_term .+= vterm
        end
        particle_list[i].force = pressure_gradient + visc_term
        #println("pressure gradient: ", particle_list[i].pressure_gradient)
    end
    
    # adj force
    accel_adj = Point3(0.0, 0.0, 0.0)
    diff_min = Point3(0.0, 0.0, 0.0)
    diff_max = Point3(0.0, 0.0, 0.0)

    map(particle_list) do particle
        accel = particle.force * mass

        speed = norm( accel )^2;
        if speed > limit^2
            accel *= limit / sqrt(speed)
        end
        
        diff_min = 2.0 * radius .- ( particle.position - pos_min ) * simscale
        diff_max = 2.0 * radius .- ( pos_max - particle.position ) * simscale
        
        # Z
        if diff_min[3] > epsiron
            normal = Point3( 0.0, 0.0, 1.0 )
            adj = extstiff * diff_min[3] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end
        if diff_max[3] > epsiron
            normal = Point3( 0.0, 0.0, -1.0 )
            adj = extstiff * diff_max[3] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end

        # X
        if diff_min[1] > epsiron
            normal = Point3( 1.0, 0.0, 0.0 )
            adj = extstiff * diff_min[1] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end
        if diff_max[1] > epsiron
            normal = Point3( -1.0, 0.0, 0.0 )
            adj = extstiff * diff_max[1] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end
        
        # Y
        if diff_min[2] > epsiron
            normal = Point3( 0.0, 1.0, 0.0 )
            adj = extstiff * diff_min[2] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end
        if diff_max[2] > epsiron
            normal = Point3( 0.0, -1.0, 0.0 )
            adj = extstiff * diff_max[2] - extdamp * dot( normal, particle.velocity )
            accel += adj * normal
        end
        
        accel += accel_g
        #
        veloc = particle.velocity
        veloc = veloc .+ accel * DT
        #
        posit = particle.position
        posit = posit .+ veloc * DT / simscale
        particle.position = posit
        particle.velocity = veloc
            #println("position, velocity: ", particle.position, particle.velocity)
    end
    positions = map(particle_list) do particle
        particle.position
    end
    positions
end

# PointsCloud settings
color = RGBA{Float32}(0.0, 1.0, 0.0, 1.0)
material = PointsMaterial(color = color, size = 0.05, vertexColors = 0)

N = 200
DT = 0.004

for i = 1:N
    println("loop: ", i)
    positions = Simulation(DT)
    positions .*= 0.1

    setobject!(vis, PointCloud(positions), material)
    trans = Translation(0., -1, -0.5)
    rot = LinearMap(RotYZ(0.2, -0.2))
    transform = compose(rot, trans)
    settransform!(vis, transform)
    sleep(0.01)
    #Base.prompt("continue?") # for capture
end