計算コードの誤り

先日の、どこが改良されたのかわからない Fortran のコードを読んでいました。
http://phreeqc.blogspot.jp/2013/10/avsmicroavsmvs.html

説明書もコメントもありませんので、数人の書いたコードを順に比較するしかありません。

最初の内は、軽微な修正で確認も早く終わりました。が、最後の方で、いくつか重要な修正がありました。変数の追加や配置、計算内容を見ると、何をしようとしているのか意図はわかります。が、どうも式がおかしい。

半日くらい考えて、たぶん間違えられたのだろうと、と結論付けかかったころ、コードが入っているフォルダの中にReadMe.txtを見つけました。お、っと思って開いてみると、ビンゴ。変更点の概略と参考文献が書かれてました。


文献を確認したが、やはりコードが違っていました。危ういですね。

といっても、その方が書かれた前の段階で、既に誤っている箇所がありました。計算部分ではなく吐き出しの箇所なので、修正は簡単です。
私の使用している post で読みやすいように、あわせて修正しておきましょう。

地震関連の調査

先日、地元の方に事業説明をしていたのですが、事前の情報が誤って伝わっており、「活断層の調査かと思っていた」と言われました。近いものはあるのですが、活断層調査ではありません。

そういえば、以前お世話になっていた方に原発の調査 = 活断層調査に声をかけられました。が、無理。おもしろそうだったんですけど。


全国で地震・耐震がブームになっているようです。

ブームで終わらせないたくないですよね。
点検や補強、数十年耐える物を作るのも良いですが、それだと一過性になりそうです。ヒトは必ず忘れます。できれば子供達への教育にお金を回したいですよね。
ベースアップした技術、というよりも意識で、さらに次の世代へつないでくれることになるでしょう。

崩壊前の地形

時間があったので、崩壊後、土砂がどの程度移動するかを再現していました（今度は LSFLOW です）。

最初は簡単に考えていました。

崩壊前後の LP を使用すれば、モデルは簡単にできてしまう、と考えていたのです。が、これが誤り。多くの場合、崩壊後の LP はありますが、崩壊前はありません。写真から復元したとしても、その精度はLPに比べて落ちます。場合によっては10m以上の標高差が生じます。今回も部分的に5～10m、最大20m近いギャップが生じました。

http://phreeqc.blogspot.jp/2012/11/5m_4.html

http://phreeqc.blogspot.jp/2012/11/5m-2.html

崩壊後のLPにしても、そのままではすべり面として使えません。冠頭部はそのままで使えるのですが、末端部は土砂で隠れています。そのため、ボーリング結果を使って面を作ってやらないといけません。こんな当たり前のことも、手を動かすまで気づきませんでした。迂闊でしたね。

こういった前処理に手間をかけることが、精度を大きく向上させるための一因になります。

案外、パラスタより時間がかかるかもしれません。

動的配列

EXCEL2010 の VBA で、TXT の並び替えをしようとしました。

あるソフトの吐き出した結果が、以下のような並び。
X1 Y1 Z11
X2 Y1 Z21
・　・　・

それを、並べなおして、新しい TXT ファイルにしたかったのです。
1 1 Z11
1 2 Z12
・　・　・

TXT ファイルを開いて、ソートして（他にもいくつか計算過程を入れて）、書き出すという流れにしようかと思ったのですが、X1, Y1 を1から始まる番号で書き直すので、Zのみの2次元配列でこなした方が早いことに気づきました。で、以下のように記述。

Dim Z(1 To NX, 1 To NY) As Long

ところが、「定数式が必要です」とエラーが出ます。
調べてみると、通常の2次元配列に2変数が使えない仕様のようです。知りませんでした。
動的配列の場合は、以下の通り書けばよいようです。

ReDim Z(1 To NX, 1 To NY) As Long


微妙です。が、忘れそうなので書き残しておきましょう。

物理演算の手抜き

3ds MAX を触っていた時に、NVIDIA の PhysX が使えるということを知りました。

MAX で PhysX の何が使えるのか、精度はどうなのか、など、知りたいことは多いのですが、なかなか詳しい解説に出会えません。が、GPU を計算に使用でき、レンダリングも可能な環境が既に用意されているのなら、使わない手はありません。少しづつ調べましょう。

先日の PBF は PhysX に取り込まれるかもしれません。
http://physxinfo.com/news/11109/introduction-to-position-based-fluids/
http://phreeqc.blogspot.jp/2013/09/blog-post.html

粒子法の SPH は既に実装済みのようです。でも、MAX でこの SPH が使えるのかは知りません。
http://sa08.idav.ucdavis.edu/CUDA_physx_fluids.Harris.pdf

ただ、引っかかっているのは「手抜き」が入っているのかどうか。以前見た本で書かれていました。
http://phreeqc.blogspot.jp/2010/11/cg.html

先日見た工学社「物理エンジンPhysX&DirectX10」にも、以下のように書かれていました。

「物理エンジン」は、与えられた時間内に処理を完了するため、必要な計算を捨てることがあります。

物理計算を”手抜き”します。

「物理エンジン」の物理学は、不完全な場合が多くあります。

でも、どのエンジンが、計算中のどこを、どの程度「手抜き」するのかが分かりません。
こういうのは、CG分野で学ばないとわからないのでしょうか？

個人的には、物理法則に忠実に問題を解いた方が better だとは思います。が、扱っている多くの問題では、所詮、パラスタで再現しているだけ = パラメーターの物理的根拠がそれらしければ良い = 同定値が真理かどうかは別問題、なので、多少手抜きしていようが再現できるのであれば良いのではないかと思います。
ま、話す場所と、相手は選ばないといけませんが。

こういうのを読むと、手を出そうかどうか、迷いますね。

微分形と積分形 その2

図書館で流体力学の本をいくつか借りてきました。

今回、役に立ったのは以下の本です。

松本洋一郎監修・山口浩樹著「道具としての流体力学」

基本的な式の展開が丁寧に掲載されています。どちらかというと、専門書というよりは一般向けかもしれません。そのため、全体的にやや物足りない感はありましたが、知識の欠落していた部分は簡単に補うことがでいました。

特に、積分形の併記が良かったですね。
質量保存則や運動量保存則より PDE を導く際、私は今まで微小領域の収支より展開する形しか知りませんでした。この本では、積分形も併記されており、2つの方法があることを知りました。これはラッキーでしたね。最近読んでいた論文が、この積分形の展開を使用していたのですが、なぜこのような数学的な展開をするのだろうと思っていたのです。これも一般的なんですね。基礎が抜けているというのは、恐ろしいことです。
微分形は微小領域（直方体）の収支から、テーラー展開経由で PDE を構成して行くのに対し、積分形は直方体にこだわらない領域の質量や運動量を体積分や面積分を使って表現し、展開して行きます（最終的には微分形で同じ式にできます）。微分形・積分形というのは、単純に式の形だけでなく、このような2種の展開の前提をも指すのかもしれません。


ナビエ-ストークスの解釈も問題ないようでした。導出も簡単です。

これで、ようやくスタートラインに立てました。

サラリーマン技術者の報酬

私の長期出張が終わったのと、次の仕事を一緒にする方が支社に来られていたのとを兼ねて、数人で飲んでいました。

いろいろ話をしていたのですが、驚いたのは報酬の話。

つい先日、私の基本給が2年間で6千円しか上がっていないことに気づきました。ここ2年間、ノルマの1.3～2倍の成果を上げ、利益を数千万出していたにも関わらずです。驚きでした。10年前と比べても、下がっていました。

これを肴に飲んでいますと、「それはまだ良い方」とのこと。　？とおもって聞いていると、皆、千円上がっただけだったり、下がったりとのことでした。世知辛い話です。

今年度の12月末までの作業量を数えると、ざっと100人。一緒に飲んでいる方との次の仕事を入れても、心身ともに余裕でこなせる量です。しかし、許可された残業＋就業時間ベースで考えると、80人程度しか費やせません。数字上は業務過多となっています。さらに先日、管理職より追加の仕事を打診されました。が、断りました。受けたとしても、サラリーマン技術者としての報酬は0円です。管理職による部下の管理もこの程度ということです。（ま、ここ2年の評価（報酬）に気づいていなければ、受けていたでしょうが。）

先日、NHK でブラック企業の話題が出ていました。

http://www.nhk.or.jp/gendai/yotei/index_yotei_3403.html

人の世です。それほどきれいな世ではありません。

F の構成

圧力項は、浸透流の動水勾配の項と似ています。

勾配が正なら逆方向に水が流れます。離れた位置の水頭が高く、動水勾配を正とすれば、逆方向に水が流れてくるので負。力も負。だから勾配にマイナスがついていると解釈すれば良いでしょう。速度勾配が大きい（正）ほど、逆向き（負）の流れが強くなるということを示している項です。

粘性項は先日の通り。プラスということは、抵抗ではなく、速度差によって微小体積に発生する力ということでしょう。同じ速度なら、ネバネバしている方が力を受けやすいというのは感覚的にOKです。
違いますね。粘性項は速度勾配の微分（速度の2階微分）なので、プラス方向の速度の差（右に流れるのを正とすれば、流速分布が右に凸な状態）が大きければ大きいほど、負になります。つまり、言葉の通り、内部の粘性抵抗ですね。

外力はプラスで問題ないでしょう。重力がZ軸の負に働けば、負の方向に速度が大きくなります。


これでナビエ-ストークスの右辺に示された3項について、特に疑問は発生しなくなったわけですが（といっても正解を見たわけではありませんが）、この3項だけで F が構成されるというのがよくわかりません。すっきりしませんが、ま、流体の多くはこれで説明できるので、そういうものだうことでしょう。

これで、SPHの前提も、今読んでいる論文の前提もイメージできました。

ラグランジュ微分のオイラー記述

∂u/∂t + (u・∇）u

ナビエ-ストークスの左辺は上記のように書かれています。これについては下記のサイトがよく整理されていると思います。
http://sky.geocities.jp/ima_ich/

Lagrange 微分の Euler 的表現
http://www2.kobe-u.ac.jp/~iwayama/teach/gfd/2012/chap2.pdf


速度のラグランジュ微分 Du/Dt （加速度) をオイラー（空間）記述しているので、移流項が出ているということでしょうねえ。この辺は、数年前に習いましたので、たぶん、あっているでしょう。

昨日の粘性項と合わせれば、移流拡散と同じ型です。スカラーかベクトルかの違いはありますが（大きな違いですが）。さらに、ナビエさんは圧力項がついているというだけのようですね。

そうすると、ナビエ-ストークスの式全体としては、加速度=力の関係になっています。右辺に1/ρがついていますが（粘性項は、動同粘性係数=粘性係数/ρとして消えていますが）、これを左辺にかけ戻せば、密度×単位体積にかかる加速度と考えられますので、ρ×V = m で ma=F、つまりニュートンの運動方程式を記述しているということでしょう。

この F として、粘性、圧力、外力を考えているということですか。どうやらここがポイントのようですね。

粘性項

ナビエ・ストークスの式で、まず引っかかったのが粘性項。

ν∆u

数式としては移流拡散の拡散項と同じ形ですし、スカラーがベクトルになっているだけです。ちなみに、2階編微分の出所は、拡散と同じで微小領域の差を取っているだけでしょう。個人的には、Dtransuのマニュアルの絵が思い出しやすいですね。

導出としては、たぶんこうだろうと、予想がつくのですが、イメージができません。拡散量はxi方向の濃度勾配に比例定数Dがかけられたもの、粘性の場合はxi方向の速度ujの勾配に比例定数がかけられたもの。うーん。どういうことでしょう。せん断方向の速度差が大きくなるほど大きくなる力が粘性ということでしょうか？

ググってみると、わかりやすい説明がありました。感謝です。
http://hb3.seikyou.ne.jp/home/E-Yama/Fluid7.pdf

運動している流体で互いに接している層の間にズレが起これば接線応力（せん断力）が現れるような流体を粘性流体と呼びました。 

ニュートンの仮定壁面よりy, y + dy の距離における速度をそれぞれu, u + du とするとy + dy にある層がy にある層をx軸方向に引っぱる力F は，この2 層の接触面積A とdu に比例し，dy に逆比例する。式で表すと比例係数をμ として

F / A=μdu/dy

これらの仮定はいろいろな実験結果より正しいことが認められています。この仮定に従う流体をニュートン流体と呼びました。
右辺F/A は，単位面積当たりの力でこれをτ で表すと

τ = =μdu/dy

となり，τ をせん断応力（Sheering Stress），比例係数μ を粘性係数（Coefficient of viscosity）と呼んでいます。また，粘性係数μ を密度ρ で割ったものを

ν=μ/ρ

で表し，ν を動粘性係数（Kinematic viscosity）と呼んでいます。

なるほど。定義なわけですね。
引っ張る速さの差が大きくなるほど大きくなるのがせん断力（抵抗力）、せん断力が現れるような流体が粘性流体、それが線形で成り立つ流体がニュートン流体、その係数が粘性係数ということでしょう。いたってシンプル。
ハチミツのように、ネバネバした液体をスプーンで早くすくうときに受ける抵抗力のようなものでしょうか？早くすくうほど、受ける抵抗が大きくなりますので。水の中で動く場合に受ける抵抗も同じですね。飛行機が飛ぶときに受ける抵抗も？要は、ネバネバして粘性係数の大きな流体や、流体中でなんらかの速度差が大きく生じる場合に取り扱うべき項と解釈しました。（流体力学の本を読まないと正解はわかりませんが）

というか、これ、昔、習っていますね。完全に忘れていますが。

ま、イメージはできたので、良しとしましょう。

微分形と積分形

式を追っていくと、「これを微分形で表すと」という表記が出てきました。

がっつり体積分の式だったのですが、いきなりそれを微分形に直すというのは意味が分かりません。過程も出ておらず、結果しか示されていませんので、変換のヒントすらありません。
そもそも、私が不勉強のため、「微分形」という用語（とその変換）があることすら知りませんでした。もっと、若いころに物理や数学を学んでおけばよかったと、改めて後悔。

で、ググってみますと（出張中なので参考書はありません）、比較的よく使われる用語のようでした。見たところ、右辺・左辺ともに体積分で統一したら、中身は一緒でしょ、その中だけ取り出せば、∇があるので微分形です、というようなことだと勝手に理解しました。ま、∇でも∂でも微分があれば微分形と呼ぶのでしょうねえ。真の定義はわかりませんが、直面している式は、その手で解決できそうです。

左辺が偏微分のついた体積分、右辺がF・nを含む面積分です。右辺をガウスの発散定理を使って体積分に上げてやると、∇・Fが生まれて、中身が一緒でしょ、はい、左辺に∂/∂tのつく微分形のできあがり。ということで、展開の結果は合いました（用語の使い方はよくわからないので、帰ったら調べないといけませんが）。

次数を上げて解きにくくしているように見えるのに、微分形の簡単な式になる、という展開の発想がすごいですね。面白い。また、物理的な意味は同じで微分形・積分形の両方で表現できるというのも（私にとっては）新鮮。

ナビエ・ストークスまでは、もう少し時間がかかりそうですが、ゆっくりいきましょう。

ガウスの発散定理のイメージ

論文を読んでいると、ナビエ・ストークスの式が出てきました。

恥ずかしながら、今までナビエ・ストークスを使ったことがなかったので、導出を見てみようと思い調べ始めました。

が、最初に出てきた質量保存式の意味が分かりません。いえ、今まで成分を書き下した式を見慣れていたので、ドット・プロダクトを使った表示の物理的意味がスッと入ってきませんでした。で、逆に成分を書き下して、体積をかけて、流速（流量）の差を取る式に直して・・・・と、イメージできるまでさかのぼって確認しました。やはり使わないと忘れますね。

以下も忘れていました。

発散 (div)   ：∇をベクトルに、形式的に内積として作用させる。∇・u
勾配 (grad) ：∇をスカラーに、形式的にスカラー積として作用させる。∇φ
ガウスの発散定理：∫∫∫v ∇・u dV = ∫∫s u・n dS

質量保存式で流入量と流出量の差を整理すると、∇・u （= ∂ui/∂xi） に行き着く。これは微小体積の発散（湧き出し量）。ガウスの発散定理の左辺は、それを（求めたい領域で）体積積分したもの。u を法線方向に直すと、微小面積当たりの流出量。それを領域の表面で面積分して湧き出し量にしたのが右辺。当然、左辺と右辺は同じ。大雑把にはこんなイメージで覚えておくと忘れにくいでしょうか？

忘れるでしょうね。若いころに叩き込んでおくべきでした。
ナビエ・ストークスは後日です。

AVS MicroAVS MVS

古い土砂移動のコードを動かしていました（LSFLOWではありません）。

数人が改良を重ねているようなのですが、どれが最新版で、どこを改良したのかコメントがありません。ただでさえコメントのないソースは、イマイチよくわからないのですが、とりあえず dimension や変数を拡張し、大きなモデルでも動くようにしました。


結果がAVSで読める形で吐き出されます。
しかし、タイムステップ毎に fld データとして吐き出されているため、アニメーションにするには一手間必要です。

まず思ったのが、AVSの利用。
ですが、他部署の方が休みのため、使用できず。

次に考えたのが、MicroAVS。
最新版がないので、試用版を落とし試しました。
補助ツールとエディターで、時系列のfldを作ってしまえばOKです。これは簡単に動画まで持っていけました。

次、MVS。
AVS の拡張が EVS, MVS なので、吐き出された fld が load EVS field　モジュールでそのまま読めるだろうと。しかも、モジュールを組んでおけば、MicroAVSより手間がかからないのではないか、と思い試すことに。
しかし、吐き出された個々の fld データや、タイムステップを記述した fld データは読めるのですが、アニメーションにできません。TCFを作っても、ステップを指定すると落ちてしまいます。当然、MicroAVS の v ファイルは読めません。ヘルプを見ても、fld を順番に読み込んで動かす方法がわかりませんでした。

うーん。ま、今日はここまで。
とりあえずはソースを動かし、アニメーションにできたので、ゴールにはたどり着けそうです。

山の粘土

離島での調査をしていました。

今回は山の調査なのですが、なぜか海成粘土のような盛土が出現。海の粘土を締固めて盛るというのは聞いたこともなく、しかし、山の粘土でもなさそう。

迷っていたので研究職に相談したところ、「分析よりも、残渣を見ろ」とのこと。

確かに。コア観察をサボってはいけません。


で、蒸発皿を現場に持っていき、コアを取って水洗い。ルーペで観察を続けましたが、貝殻片や黄鉄鉱は入っていません。結論から言うと、海のものではなく、淡水のもの。

しかし、どこから、どのように堆積した粘土を持ってきたのでしょうか？結局、そこまではわかりませんでした。


常日頃、1m下の地盤や岩盤に振り回されているレベルです。まだまだです。