趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

ライセンスの話

例のSection Quasi 2Dモデルについて、Githubで公開しようかなと考えているのですが、そろそろライセンスのことも気にしたほうが良いかもしれないと思ってます。

どうしましょうか。MITかGPLかと思いますが、GPLにするほど大したものじゃないですが悩ましいです。

作者はこのHPにします。ミスがあればその都度gitで直していこうかと思います。

また、jupyterの部分はpythonのコードに変更しようと思います。GUIは作らないです。

流砂の話とか

  • fortranpythonですが、calSectionProfileとinitializeを呼んでます。いろいろ試していると、fortranの共有部分でallocateしたものが、ずっと残ってしまうので、計算が終わるごとにdeallocateしたほうが良さそうです。また、まとめて報告します。

  • XeonGPUとか使いたいのですが、まだまだ計算の力が追いついてないです。すごいコードができから回したいですね。

  • 今年最後になると思いますが、ゆっくりだけどそれなりコードを書いた一年でした。徐々にストックも溜まってきてうれしいです。来年も頑張ります。

断面準二次元モデルの修正方針

  • 自分の中での反省ですが、自分の中では粗度係数はパラメータじゃないと言ってきたのですが、このモデルでは相当粗度を同化するとか言ってるので矛盾していることに気付きました。
  • そこで、相当粗度は、例えば、黒木・岸のΦ-τ*の関係等から断面平均的に小規模河床波に応じた抵抗を設定し、それより流速場を計算し、観測結果との比較により、duneっぽいなとかを議論できるモデルに変更しようと思います。

流砂の話

  • 最近感じるのですが、流砂はDEMじゃ無理だな思います。というのも、河川流では粘性底層の流れでさえも未だに解けていないため、掃流砂や浮遊砂のように明らかに流れ場に影響を与えるようなものに対して、その間隙流を解析することは不可能だと思うからです。
  • 仮に河床にNon-Slip条件を与えて解けるモデルがあれば、粒子の面にもNonslip条件を与えて(IB法とか面白そう)、流れ場とのインタラクションが評価できるので可能と思いますが、流体力学の分野でもまだまだなので、ちょっと遠いかなと思います。
  • 逆にwashloadのような流れ場に影響与えないものはDEMでも可能かと思います。
  • 当面は、掃流砂・浮遊砂は連続体しかないと思ってます。昔流行っていた二層モデルなんかが良いような気がします。
  • 先日ボスと話していたら、流砂研究が人気のない要因は、意味の分からない経験則ばかりだからだと言っていて、確かに流砂の運動方程式がないのはまずいなと思いました。微分方程式がないので、例えば計算で精度を上げようとしても無理ですし、楽しみもないです。流砂は、河道管理の最重要課題なので、何とかしたいものですね。

平面二次元モデル

  • 数ヶ月前の話ですが、例のコロケート格子で流速0の場合に凹凸が消せない件ですが、何かアイデアは浮かびましたでしょうか。私は、どう考えても分からないので非物理的なフィルタで消すことしか思い浮かびません。次回にでも相談させて下さい。

FortranとかPythonとか流砂とか

Fortran + Python

Fortranとの連携

ソースコードだけ見ました。力作ですね。

fortranのcalSectionProfileを呼び出してるということで良いんでしょうか?

Pythonを使ってしまうと、行列の演算は他の言語だとキツくなってきますね。 fortranもライブラリを使えば良いのでしょうが、ちょっとしたものだとどうしてもPythonですね。

楽なのでPythonで書いてしまうのですが、2次元以上だと計算時間が厳しくなって、 結局Fortranで書き換えるのですが、辛いです。

あともはや最近はCを書けなくなってます。ヘッダが書きたくないです。

ソースはWindowsで試してみますね。

matplotによる可視化

可視化が使えるのが良いですね。グラフも綺麗ですし。

仕事では、gnuplotが多いですが、matplotの方がカッコイイです。

会社の人間はあまり新しいプログラミング言語だったり、新しいツールに興味に持つのが難しいようですね。 仕事なので仕方ない部分もありますが。

フリーのコンパイラ

先月PGIのコンパイラにCommunity Editionが出た模様です。

https://www.softek.co.jp/SPG/Pgi/product.html

基本はフリーで使えますが、最新版のみのライセンスのようです。

また、LinuxMacのみでWindowsは対応未定のようです。 PGIはGPUとの連携が強いので、使ってみたいです。

Macでは入れてみましたが、特に問題なく動きました。 動作チェックしかしていないので、実行速度については見ていませんし、 家にGPUマシンがないので、ちゃんと使えてはいませんが。。。

基礎水理シンポジウム

  • 村上さんの発表は現場からの話ですが、土石流と河床変動をシームレスに解析するのは、面白いですね。本来別のものではないはずなので、江頭先生や伊藤先生の論文のように現場でもそういう解析やるのは大事だと思います。

  • 浮遊砂を含めて、今のフレームワークを取っ払ってゼロから考えないとダメなのかなと思ってます。あとはこれらの研究が実用上いかに重要かをどう説明するかですね。

基礎水理シンポジウムとか

昨日所用で北海道に行ってました。寒すぎて生きている心地がしなかったです。

f:id:SedimentHydraulics:20161217185120j:plain

基礎水理シンポジウム

先日行われた基礎水理シンポジウムで気になった内容をまとめておこうと思っていましたが、すっかり忘れておりました。

  • テーマは、浮遊砂・ウォッシュロードの力学についてです。発表者は、江頭先生、村上さん(CTI)、内田さん(国総研)、権田先生、堀田先生、関根先生です。部会長が里深先生ということもあり、砂防よりの人選です。
  • 後日PPTが公開されると思うので、詳細はそこで確認いただければと思いますが、全体的にかなり真面目に浮遊砂問題にアプローチしているなと言う印象です。河川も見習うべきだと思いました。
  • 皆さんの話を聞いていると、平衡状態では、浮遊砂の輸送過程の力学は理解できているのかなと感じました。ただし、浮遊砂の境界条件については全然駄目っぽいです。特に、浮遊砂層の底面位置(いわゆる基準面濃度の位置)とそこからの巻き上げ量については、力学が無いような状況です。
  • 基準面濃度については、us/w0で整理するのは無理じゃないかという話がでてました。確かに計測値の誤差が大き過ぎて評価が難しいです。
  • 底面の位置は、掃流砂層の上面になるので、安定・不安定解析みたいな方法で決めれないかなと言う話があり、面白そうだな思いました。なんとなくですが、乱流みたいに掃流砂層、遷移層、浮遊砂層みないな感じなるのかなと思います。
  • また、us/w0が1付近は、少し複雑な力学になるようです。河川の世界では、主になる領域なので重要かと思います。
  • 高濃度流れでは、沈降速度が粒子単独のものよりは、小さくなるという実験結果もありました。ウォッシュロードが卓越するような場では考える必要があるのかなと思います。

まだまだ面白くなりそうです。

Section Quasi2D model in natural river

できました。

とりあえず、結果だけ。いい感じです。

f:id:SedimentHydraulics:20161210163514p:plain

ptyhon+fortranです。ちょっと力作なのでコードを見てほしいです。

DBに置いておきます。

python+fortran:2003+OOP

本業が忙しくて全然だめです。

python+fortran:2003+OOP

以前書いた話の続きですが、pythonからfortranのdllを使うときもやっぱりOOPで書きたいので、チャレンジしてみました。

結果から言うと、main文に共有moduleを設け、その中で共有のclassを設定することにより雰囲気はでました。 C系だとpythonからobjectを渡せるのでこんなことをする必要はないと思いますが、fortranだとこれが限界です。

あと忘れないように書いておくと、fortranのfunctionは使えません。

!file name : test.f90
!class
module classModel
    implicit none
    private    
    type, public :: model
        DOUBLE PRECISION :: x
    contains
        procedure, public :: init
    end type
contains

subroutine init(self, x)
    class(model) :: self
    DOUBLE PRECISION, INTENT(IN) :: x
    self%x = x
end subroutine

end module classModel

!main
module pubMod
    use classModel
    CLASS(model), ALLOCATABLE :: m
end module

subroutine init(x)
    !DEC$ ATTRIBUTES DLLEXPORT :: init
    use pubMod
    DOUBLE PRECISION, INTENT(IN) :: x
    allocate(m)
    call m%init(x)
end subroutine

subroutine plus(y, ans)
    !DEC$ ATTRIBUTES DLLEXPORT :: plus
    use pubMod
    DOUBLE PRECISION, INTENT(IN) :: y
    DOUBLE PRECISION, INTENT(OUT) :: ans
    ans = m%x + y
end subroutine
@echo off
gfortran -std=f2003 -shared -o test.dll test.f90
import numpy as np
from ctypes import *

f = np.ctypeslib.load_library("test.dll",".")

f.init_.argtypes = [POINTER(c_double)]
f.plus_.argtypes = [POINTER(c_double),POINTER(c_double)]

x,y = 2.0 ,1.1
x = c_double(x)
y = c_double(y)
res = c_double()

f.init_(byref(x))
f.plus_(byref(y), byref(res))
print(res.value) # 3.1

del f

DBに置いておくので回してみて下さい。

これを使ってsection2Dモデルの拡張版を書いているので近々公開します。

WindowsFORTRAN環境構築について訂正

先日gfortranの4.8が良いって言いましたが、6.0代も安定してますのでこっちのほうが良さそうです。

64bit の Windows10 上でフリーの fortran コンパイラを導入して、簡単なプログラムを作成する - あらきけいすけの雑記帳

exe一発です。

windowsfortranコンパイル用のmakefileって....

最近知ったのですが、inteのnmakeとgmakeって結構違いますね。どちらでもmake一発でという思いでmakeファイルを書いていたので、残念です。

書く気をなくしました。

不等流計算と不定流計算

結構忙しくて全然書けておりません。すいません。

不等流計算と不定流計算について

メールでも送った話ですが結構面白かったので、書いておきます。

不等流計算と不定流計算の定常解の計算値に結構差がでるなと気付き、色々考えてみました。

計算結果がフィールドを対象としているため、結果は出せないのが辛いところです。

物理的な根拠としては、Frが1以下の場合は下流の影響のみで評価できるというのは矩形で、dB/dxが漸変的な流れで導出される考え方であり、当然上下流の影響は受けるとだろうということだと思います。

また、一般断面の場合、Frの定義にもdA/dBを考慮する必要があると思います。

後日式関係も整理しておきます。

この辺のからみで高橋保先生の面白い実験があったので、再現計算してみようと思います。

いっしょにやりませんか。

一般断面のHLLCを書いてほしいなあ。。。

FORTRAN環境ってあります?

たしか以前Mac OSのアップデートでgfortranに問題あるって言ってましたけど、環境はどうでしょうか?

上記のコードはF2003で書いているのでできれば環境を構築してほしいなと思ってます。gfortran4.8シリーズが安定してるのでおすすめです。5.0シリーズは全然使えないです。

関係ないですがVScode使い始めました。いい出来ですね。まだ慣れてないですが、言語に よってはデバッグもできるようだし、IDE的にもなるのかなと思ってます。