#!/usr/bin/env ruby
# -*- coding: utf-8 -*-
# 
# 1 次元定常移流方程式
# \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0
#
# 時間方向に 4 次のルンゲクッタ法, 空間方向に 4 次の中央差分を取った場合.
# 実習用ではなく, 講師側の確認用
# 
# = 履歴 (編集した際は上に追記すること)
# * 2012/06/12 Satoshi NODA  作成
# 

require "numru/ggraph"
include NumRu
include GGraph

### パラメータ設定
nt     = 30000      # 全ステップ数
nt_out = 1000       # 出力間隔

nx = 128           # 格子点数
l  = 2 * Math::PI  # 系の大きさ
c  = 1             # 移流速度
dx = l / nx        # 格子の間隔
dt = 0.001         # 時間ステップ

### 座標の設定
x = dx * NArray.float(nx).indgen

### 初期値の設定
u0 = NArray.float(nx).fill(0.0)
n1 = nx / 8
n2 = nx / 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
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 draw(u, x, t)
  # gphys 形式へ変換 (NetCDF 形式で出力するときに便利なように)
  gphys = na2gphys(u, x)
  # 折れ線グラフとして可視化
  line(gphys, true, 'max'=>1.7, 'min'=>-0.7, 'title'=>"t="+t.to_s)
end

def dudx(u, dx)
#  # du / dx (2 次の中心差分)
#  dudx = (u.sft(1) - u.sft(-1)) / (2*dx)
  # du / dx (4 次の中心差分)
  dudx = (- u.sft(2) +8*u.sft(1) -8*u.sft(-1) + u.sft(-2)) / (12*dx)
end

########################################
# メイン
########################################

### ウィンドウを開く
DCL.gropn(4)

### 初期値の代入と描画
u = u0
draw(u, x, 0.0)

### メインループ開始
(1..nt).each{|it|
  # 時間発展を解く (4 次のルンゲクッタ法)
  k1 = dudx(u, dx)
  k2 = dudx(u+0.5*dt*k1, dx)
  k3 = dudx(u+0.5*dt*k2, dx)
  k4 = dudx(u+dt*k3, dx)
  dudt = -c * (k1 + 2*k2 + 2*k3 + k4) / 6.0
  u = u + dt * dudt

  # 可視化
  if it % nt_out == 0
    t = it * dt
    draw(u, x, t)
  end
### メインループ終了
}

### ウィンドウを閉じる
DCL.grcls