Fantom Refined Manual

\(\require{physics}\)

Particle Trackingオブジェクトの使い方

New

各パーティクルトラッキングが出力するVTKファイルのファイル名（ヘッダ部分）を指定できるように変更した（今までは固定であったため，ファイルが上書きされてしまっていた）．基本的にオブジェクトの生成の際の引数の最後に名前を指定する．

particle_tracking = particle.LiveParticleTracking(water_body, "your_header_name")
particle_tracking = particle.PlaybackParticleTracking(water_body, files, input_interval*dt, "your_header_name")
particle_tracking = particle.PlaybackParticleTracking(water_body, files, input_interval*dt, string_vector({"salinity"}), "your_header_name")

概要

数値モデルを用いて水域内の物質移動を解析する手法として，トレーサー（色素のようなもの）による解析とパーティクル（粒子）による解析が考えられる．それぞれの解析方法には得手不得手がある．トレーサー解析の場合，トレーサーの移流・変形を十分にグリッドで解像できない場合，数値拡散が支配的になってしまう可能性がある．例えば，大領域における物質の移動先の特定を主目的とした解析，例えば，粒状性物質（卵など）の輸送解析には不向きである．一方，パーティクルの解析では，個々の粒子の位置は具体的に追跡できるが，性質上，連続的な表現ができない．汚染物質の拡散などの濃度解析を考える場合に，多数の粒子で解析することで濃度も評価できるかもしれないが，粒子の個数の増大に伴う負荷増大が問題になる．

パーティクル利用のメリットとして，パーティクルの運動が流れに影響を与えない場合（パッシブ・トレーサーの場合）は，計算済みの流れデータを用いて粒子の移動を計算できることが挙げられる．Fantom Refinedでは，流れ解析と同時に粒子の追跡計算(LiveParticleTracking)を行うこと，事前に計算した流動の結果(datデータ)を読み込みながら，流体計算をせずに粒子追跡を行う方法（PlaybackParticleTracking）の両方を利用することができる．追跡する際の流動場情報の時間間隔が異なれば，両者の結果は完全に同じにはならないが，ファイルで保存されている流動場が流れの時間変化を十分に表現していれば，粒子の集団としての挙動で大きな差が出ることはない．

大まかな手順は，以下の通りである．

ParticleTrackingオブジェクトを作成する．
ParticleTrackingオブジェクトに粒子を任意の数だけ，任意の位置に追加する（動的に流動計算中に追加しても良い）．
時間進行ループの中でParticleTrackingオブジェクトをアップデートして粒子の移動を行う．
適当なタイミングでvtk形式で出力する（もしくは，Luaの簡単なコードを書いてファイルに出力する）．
粒子の沈降速度などを考慮する場合は，各パーティクルに沈降速度を指定する．

オブジェクトの種類と注意

LiveParticleTracking：流れと同時に粒子追跡を追跡する．
PlaybackParticleTracking：既存の流動計算結果を用いて粒子追跡を追跡する．
Particleを用いる場合は必ずLuaファイルの先頭でrequire("particle")を記述する（ライブラリーを読み込む）．

LiveParticleTrackingの利用方法

ParticleTrackingオブジェクトはWaterBodyオブジェクトを引数にするので，WaterBodyを作成した後であれば問題ないが，時間ループの直前あたりで生成するのがわかりやすいと思われる．

water_body        = fantom.WaterBody(dx, dy, dz, ....省略)
particle_tracking = particle.LiveParticleTracking(water_body)
particle_tracking:add_particle(x0,y0,z0) -- ３次元的に動くパーティクルを一つx0,y0,z0の位置に追加
particle_tracking:add_floating_particle(x0,y0) -- 水面上を動くパーティクルを一つx0,y0の位置に追加

上の例では１つ一つ粒子を追加している．つまり，同じコマンドを繰り返し使うことでメモリーの許す限りのパーティクルを追加することができる．そのほか，位置(x,y,z)座標を別途ファイルに記述し，読み込んで大量の粒子を生成することもできる（オプション）．

particle_loc_file = io.open("location_data.txt", "r")
for line in particle_loc_file:lines() do 
  data = split(line, "\t") -- split関数は別途定義
  particle_tracking:add_particle(data[1],data[2],data[3])
end

また，コンテナに満遍なく等間隔で粒子を撒く場合には以下の関数を利用すると良い．以下の関数を呼ぶ場合は，コンテナの(i,j)と粒子を配置する間隔(x_interval, y_interval, z_interval)を指定する（オプション）．

function container_particle_seeding(water_body, particle_tracking, i, j, x_interval, y_interval, z_interval)
  num_particles = 0
  west_end = water_body:dx()*i
  south_end = water_body:dy()*j
  nx = water_body:dx()/x_interval
  ny = water_body:dy()/y_interval
  nz = water_body:thickness()/z_interval  
  for i=0, nx do
    for j=0, ny do
      for k=0, nz do
        if particle_tracking:add_particle_bool(west_end+i*x_interval, south_end+j*y_interval, k*z_interval)==true then
          num_particles = num_particles+1
        end
      end
    end
  end
  return num_particles
end

流体計算を行う時間進行のループ内では以下のように粒子の移動及びファイルの出力を行う．

for counter = 0, total_count do
  time = dt*counter
  if counter % save_interval == 0 then
    particle_tracking:write(counter*dt, counter)
  end	
  ・・・
  particle_tracking:move_particles(dt)
end

PlaybackParticleTrackingの利用方法

プレイバックの場合は，まず過去のデータを計算するために利用したデータセット（地形など）を用意する必要がある．そして粒子追跡に利用する流動データファイルのファイル名リストをなんらかの形で用意する（直接Luaファイルに書き込んでも構わない）．また，プレイバックの場合は粒子の追跡時間間隔dtの他に，読み込みデータの時間間隔を指定する必要があることに注意する．以下にプレイバック粒子追跡を行うプログラム例を示す．

require("fantom")   --このライブラリーを読み込むこと
require("particle") --このライブラリーを読み込むこと

function string_vector(tbl) -- 文字列のテーブルをvectorへ変換
  local tmp = fantom.VectorString()
  for i,v in ipairs(tbl) do tmp:push_back(v) end
  return tmp
end

-- basic parameters
case_name              = "pb_example"
dt                     = 3.0                -- 粒子追跡時間間隔
save_interval          = 30*60.0/dt         -- 保存時間間隔
input_interval         = 5*60/dt            -- 入力データ間隔
dx                     = 20.0*64            -- 流動計算に用いたものと同じ
dy                     = 20.0*64            -- 流動計算に用いたものと同じ
total_count            = (2*24+2)*3600.0/dt -- 追跡ステップ数
limit_depth            = 0.05               -- 流動計算に用いたものと同じ

dz         = fantom.single_line_data("./input/dz_data.inp")
water_body = fantom.WaterBody(dx, dy, dz, "bottom_data.inp", "xrfm.txt", "yrfm.txt", limit_depth)
files      = string_vector({"data01.dat", "data02.dat"}) -- 必要なデータファイルを全て記入（非常にたくさんの場合はfor-loopで作成）

-- create playbackParticleTracking object
particle_tracking = particle.PlaybackParticleTracking(water_body, files, input_interval*dt, string_vector({"salinity"})) --粒子位置の塩分を求めるためにsalinityを追加

-- add particles
particle_tracking:add_particle(0.0, 0.0, 0.0) -- 例えば，(0,0,0)に１個のパーティクル追加

-- set settling velocity (この変更を時間進行loop内で行うこともできる)
for n=0, particle_tracking:particles():size()-1 do
  particle_tracking:particles()[n]:set_additional_velocity(0,0, 0.0, 0.0)
end

-- time marching
for counter = 0, total_count do
  print("counter = ", counter)
  if counter % save_interval == 0 then
    particle_tracking:write(counter*dt, counter) -- vtkファイルを出力
  end
  particle_tracking:move_particles(dt)
  collectgarbage()
end

addditional_velocityを塩分やその他のスカラー量に基づいて変化させることによって，卵の移動や魚の挙動をある程度模擬することができるようになる．

最後に粒子位置をカウンターの番号ごとにファイルに保存するLua関数を示す．

function save_location(particles, counter, particle_tracking)
  local f = io.open("./result_particle/location_"..string.format( "%06d", counter )..".txt", "w")
  for i=0, particles:size()-1 do
    f:write(tostring(particles[i]:loc().x).."\t"..tostring(particles[i]:loc().y).."\t"..tostring(particles[i]:loc().z).."\t"..tostring(particles[i]:current_column():surface_level().present-particles[i]:loc().z).."\t"..tostring(particles[i]:current_column():surface_level().present-particles[i]:current_column():bottom_level().present).."\t"..tostring(particle_tracking:get_particle_scalar(i,"salinity")).."\n")
  end
  f:close()
end

利用するためには保存するタイミングで以下のように上記関数を呼ぶ．

save_location(particle_tracking:particles(), counter, particle_tracking)

サンプルは以下の通り

require("fantom")
  require("particle")
  
  -- utility function
  function string_vector(tbl);
    local tmp = fantom.VectorString()
    for i,v in ipairs(tbl) do tmp:push_back(v) end
    return tmp
  end
  
  function save_location(particles, counter, particle_tracking)
    local f = io.open("./result_particle/location_"..string.format( "%06d", counter )..".txt", "w")
    for i=0, particles:size()-1 do
      f:write(tostring(particles[i]:loc().x).."\t"..tostring(particles[i]:loc().y).."\t"..tostring(particles[i]:loc().z).."\n")
    end
    f:close()
  end
  
  -- basic parameters
  case_name              = "./result/ogouchi(casex01)"
  dx                     = 6.25*8
  dy                     = 6.25*8
  dt                     = 20
  total_count            = 840
  save_interval          = 600/dt	-- every 10 mins
  input_interval         = 600/dt
  limit_depth            = 0.05
  
  dz = fantom.single_line_data("./input/dz.inp")
  water_body   = fantom.WaterBody(dx, dy, dz, "./input/geometry_50m.inp", "./input/container_50m.inp", "./input/container_50m.inp", limit_depth)
  particle_start      = 0
  input_time_interval = input_interval*dt
  
  -- file names
  data_names = {}
  for i = 0, total_count/input_interval do
      table.insert(data_names, "./result/ogouchi(casex01)-"..string.format( "%06d", math.tointeger(particle_start+i*input_interval))..".dat")
  end
  files = string_vector(data_names)
  --- particle tracking object
  particle_tracking = particle.PlaybackParticleTracking(water_body, files, input_time_interval) 
  
  --- add particles
  dist = 20 --m
  for i=0, 500 do
    for j=0, 500 do
      particle_tracking:add_floating_particle(i*dist, j*dist)
    end
  end
  
  -- time marching
  for counter = 0, total_count do
    print("counter = ", counter)
    if counter % save_interval == 0 then
      particle_tracking:write(counter*dt, counter) -- vtk
      save_location(particle_tracking:particles(), counter, particle_tracking) -- txt
    end
    particle_tracking:move_particles(dt)
    collectgarbage()
  end

また，出力結果（テキスト）の可視化スクリプトは以下の通り

import glob
  import matplotlib.pyplot as plt
  import numpy as np
  
  files = sorted(glob.glob("./result_particle/*.txt"))
  
  for file in files:
    print("reading ... ",file)
    data = np.array([list(map(float, line.rstrip().split())) for line in open(file).readlines()]).transpose()
    plt.axes().set_aspect('equal')
    plt.scatter(data[0], data[1], alpha=1.0, marker='.', c='k', s=2, edgecolors='none')
    plt.xlabel("x (m)")
    plt.ylabel("y (m)")
    plt.title("particle tracking")
    plt.gca().set_xlim(0,6500)
    plt.gca().set_ylim(0,3500)
    plt.savefig(file[:-4]+".png", dpi=300)
    plt.clf()