[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[dennou-ruby:001815] Re: GPhys::EP_flux
堀之内です。
＃ 石岡さんdennou-ruby ML 入ってましたっけ？
塚原君は作業で頑張ってますので、私が代りにメールします。
作戦変更レポート。
● SHTLIB (in DCL/MATH2にある球面調和関数展開ライブラリー) について
等間隔データが対象ということで、等間隔化して提供されている客観解
析データの処理にはもってこいということで使う話でした。RubyDCLに
含まれてますので、すぐ使えるのも魅力。ただ、RubyDCLのSHTLIBイン
ターフェースにはちょっとバグがあるのを塚原君が発見し、直してくれ
ました。次期リリースには取り込まれます。（もしすぐ使いたいという
声があればすぐにリリースしますが、今のところ SHTLIB 以外のアップ
デートはありませんので。）
[[ここで石岡さんに質問です。初期化ルーチン shtint では、
JM >= (MM+1)/2, IM >= MM+1 でなければならない
と書いてありますが、やってみると IM, JM がそれよりある程度小さく
てもうまく行くように見えます。ただし小さすぎると実際に変換→逆変
換で戻ってきませんし、下手すると segmentation fault がでます。実
際のところ、どうなんでしょう。特に、MM+1 の 「+1」は必要ですか。
とすると、例えば東西に10分割した場合に、波数が 5 まででなく 4 ま
でになってしまいます。ちなみに他のルーチンのドキュメントに少々間
違いを発見しました。そのうち塚原君から報告があると思います]]
(石岡さんへの質問おわり)
さて、SHTLIBはどうやら三角切断専用なようです。東西波数を南北より
も小さくとれないので、EP flux で出てくるような帯状平均物理量、す
なわち波数ゼロ成分だけに関する微分計算に使うには不向です。むりや
り東西方向に同じものをコピーして（つまり ary(ny) --> any(nx,ny)
のように広げちゃう）使うというのも考えたのですが、あまりに無駄が
多い(∵ nx=~2*ny)ので、やめようということになりました。
あと、下手な関数を与えるとギブス現象がでます。参考まで。ちゃんと
考察してませんが、極での特異性があるべき状態からはずれると起るよ
うに思います。よって、球面geometry的にリーズナブルな関数を与え、
かつ球面geometry的にリーズナブルな計算をしてる限りは問題にならな
いような気がします。ということで話はつきそうなもんですが、一応以
下に具体例をあげます。ちゃんと考えないとトラブるという実感が共有
して頂けると思います。
テストその１： g(λ,φ) = cosφ (帯状一様)
shts2g で isw=1 とすると ∂g/∂φ が計算できますが、極で１になる
べきところがゼロになり、そのぶん相当シビアなギブスの角がでました。
一方、shtg2s で isw=1 とすると(cosφ)^-1∂(g cosφ)/∂φ が計算
できますが、この場合は問題ありません。実際、EP flux 計算で帯状平
均量に対して適用するのは、問題のなかったほうの後者です。
テストその２： g(λ,φ) = sinλ cosφ
この場合は ∂g/∂φ はうまく行き、(cosφ)^-1∂(g cosφ)/∂φ は
激しいギブスの角がでました。
● で、微分はどうするの
球面というgeometryに特化したモジュールを作らないことになったので、
汎用の微分用モジュールをつくる。モジュール名は NumRu::Deriv か
NumRu::Derivative（対象はNArray）。また、これを GPhys に組み込ん
だものもつくる。
とりあえず、中央差分（dx_i = x_i+1 - x_i-1 タイプ）をつくる。
境界条件は取り替え可。デフォルトは dx_0 = x_1 - x_0,
dx_N = x_N - x_N-1。実装は境界条件にあわせて配列を拡張したのち、
x[2..-1] - x[0..-3] をとる。境界条件にあわせた拡張をメソッド化す
ることで、境界条件変更を容易にする。ちなみに今回の目的に照らすと、
極では上記のデフォルトのが妥当と考えられるので、とりあえずこの１
種類のみをつくる。（滞在時間が限られるので、基本方針は、「汎用化
のことを考えて設計しつつ、実装は必要最低限で突き進む。」）
> =begin
> = GPhys 用 EP-flux 計算モジュール作成作戦メモ
>
> このメモは, 京大生存圏の堀之内さん部屋で行われた EP-flux モジュール作成に関する
> 堀之内さんと塚原の打ち合わせメモである.
>
> == 履歴
> 2004/08/03(火) 塚原大輔 -- 新規作成.
>
>
> == 作成するモジュール
> 名称は仮である.
>
> * ((<module NumRu::Spherical>)) / ((<module NumRu::GPhys::Spherical>))
> * 球面座標での微分演算を行う関数群. NArray, GPhys の双方に対して作成.
> * ((<module NumRu::Cartesian>)) / ((<module NumRu::GPhys::Cartesian>))
> * 直交直線座標での微分演算を行う関数群. NArray, GPhys の双方に対して作成.
> * ((<module NumRu::GPhys::EP_flux>))
> * EP_flux を算出する関数. GPhys に対して作成.
>
> === 微分演算 モジュール
>
> * NumRu::Spherical
> * 球面座標での微分演算モジュール.
> * xderiv, yderiv, rot, div, xdiv, ydiv, grad, xgrad, ygrad, lapla, ilapla, jacobian, .., ydiv_m0, ygrad_m0, yderiv_m0, etc.
> * 等間隔グリッドデータを前提とする
> * その他のグリッド(特にガウス緯度)についてはおいおい対処.
> * ISPACK を導入しやすいよう意識して枠組みを作成
> * スペクトル変換, 微分については DCL の SHTLIB を利用
> * ((<URL|http://www.gfd-dennou.org/arch/ruby/products/ruby-dcl/ruby-dcl-doc/math2/node24.html>))
> * NumRu::Cartesian
> * デカルト座標での微分演算モジュール.
> * スタッガードに対する関数はとりあえず作らない.
> * とりあえず差分法
> * cderiv(中央差分), fderiv(前方差分)..etc.
>
> * NumRu::GPhys::Spherical
> * NumRu::Spherical を利用
> * データの定義領域のチェック(緯度は何度から何度までか, 緯度経度なのか緯度のみか, とか)
> * メソッド化して, 下請けに応じて取り替えられるようにする
> * ISPACK と SHTLIB でデータの定義領域が(ちょっと)異なる
> * 軸の単位をチェック(角度かどうか)して, 判定
> * 北から南にデータが入っている場合(NCEP はそう)は, 逆向きにしてあげる
> * 将来的には, (可能な限り)領域を適合させられるようにする.
> * 単位, 地球半径を考慮
> * NumRu::GPhys::Cartesian
> * NumRu::Cartesian を利用
> * 座標変数の自動更新
> * 例 : 半整数グリッドと整数グリッド
> * 単位を考慮
>
> === EP_flux モジュール(NumRu::GPhys::EP_flux)
> * 基本は近似を行わない, いわゆるフルセットの EP_flux
> * 近似を導入しようとも, 微分演算モジュールは作らねばならないため, コストは一緒.
> * オプションで近似形を用いるかどうか切り替えられるようにする.
> * 入力変数
> * 東西風, 南北風, 鉛直風(omega or w), 温度(temp or theta)
> * 鉛直座標は圧力座標(P or logP)
> * 出力変数
> * Fy, Fz, uv_bar, ...
> * divF 計算メソッドを作成
> * Fy, Fz を引数にとる
>
> == スケジュール(目標)
>
> (1) NArray 版微分演算モジュールの作成 -- 8/4(水)
> (2) GPhys 版微分演算モジュールの作成 -- 8/5(木)
> (3) EP_flux 算出モジュールの作成 -- 8/6(金), 8/7(土)
> (4) チュートリアルの作成 -- 8/8(日), 8/8(月)
>
> === 留意点
> * コーディングとテストは同時進行
> * テストの方針
> * 解析的解のデータと出力を比較する
> * テストがある程度通るようになったら, ドキュメントを作成(日本語 and RD)
>
> =end