#!/usr/bin/env ruby # -*- coding: utf-8 -*- # # 1 次元定常移流方程式 # \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 # # 時間方向に 1 次の前進差分, 空間方向に 2 次の中央差分を取った場合. # # = 履歴 (編集した際は上に追記すること) # * 2012/07/25 納多哲史 解析解を一緒に表示するよう変更 # * 2012/06/12 納多哲史 作成 # require "numru/ggraph" include NumRu include GGraph ### パラメータ設定 nt = 9000 # 全ステップ数 nt_out = 300 # 出力間隔 m = 128 # 格子点数 l = 2 * Math::PI # 系の大きさ c = ARGV[0].to_i # 移流速度 ARGVはコマンドラインの引数の文字列の配列 to_iは文字列を数に変換する。 dx = l / m # 格子の間隔 dt = 0.001 # 時間ステップ ### 座標の設定 x = dx * NArray.float(m).indgen ### 初期値の設定 u0 = NArray.float(m).fill(0.0) n1 = m / 8 n2 = m / 4 u0[n1..n2] = 1.0 ######################################## # メソッド定義 ######################################## class NArray ### 周期境界条件で i (整数) 個シフトさせた配列を返す def sft(i) na = NArray.float(self.size) na[0..(-i-1)] = self[i..-1] na[-i..-1] = self[0..(i-1)] na end ### 周期境界条件で xsft (浮動小数点) 個シフトさせた配列を返す. ### 解析解の計算用. def sft_float(xsft, l, m, dx) na = NArray.sfloat(self.size) sft = (xsft % l) / dx dm = sft.divmod(1.0) isft = dm[0] sft_mod = dm[1] na = (1.0 - sft_mod) * self.sft(isft) + sft_mod * self.sft((isft+1) % m) na end end ### 1 次元の NArray からおまかせで GPhys を作成 def na2gphys(nary, x) # option axis_opt = {'long_name' => 'X-axis', 'units' => '1', 'name' => 'x'} data_opt = {'long_name' => 'u', 'units' => '1', 'name' => 'u'} # 軸情報の設定 va_axis_x = VArray.new( x, {"long_name"=>axis_opt['long_name'], "units"=>axis_opt['units']}, axis_opt['name'] ) axis_x = Axis.new.set_pos(va_axis_x) # データとその付属情報の設定 va_data = VArray.new( nary, {"long_name"=>data_opt['long_name'], "units"=>data_opt['units']}, data_opt['name'] ) # GPhys オブジェクトの作成 gphys = GPhys.new( Grid.new(axis_x), va_data ) end ### 解析解の計算 def u_ref(u0, xsft, l, m, dx) u_ref = u0.sft_float(xsft, l, m, dx) u_ref end ### 描画 def draw(u, x, t, u0, l, m, dx, c) # gphys 形式へ変換 (NetCDF 形式で出力するときに便利なように) gphys = na2gphys(u, x) gphys_ref = na2gphys(u_ref(u0, -c*t, l, m, dx), x) # 折れ線グラフとして可視化 # 解析解の描画 line(gphys_ref, true, 'max'=>2.0, 'min'=>-1.0, 'index'=>15, 'title'=>"t="+t.to_s) # 数値計算によって得られた値の描画 line(gphys, false, 'max'=>2.0, 'min'=>-1.0, 'index'=>25) end ######################################## # メイン ######################################## ### ウィンドウを開く DCL.gropn(4) ### 初期値の代入と描画 u = u0 draw(u, x, 0.0, u0, l, m, dx, c) ### メインループ開始 (1..nt).each{|it| # du / dx (2 次の中心差分) if (c < 0) then dudx = (u.sft(1) - u.sft(0)) / dx else dudx = (u.sft(0) - u.sft(-1)) / dx end dudt = -c * dudx # 時間発展を解く (一次の前進差分) u = u + dt * dudt # 可視化 if it % nt_out == 0 t = it * dt draw(u, x, t, u0, l, m, dx, c) end ### メインループ終了 } ### ウィンドウを閉じる DCL.grcls