趣味で計算流砂水理

趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

数値計算とか河川工学とかプログラミングのことを書いています

MENU

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マシンがないので、ちゃんと使えてはいませんが。。。

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的にもなるのかなと思ってます。

ニ次元で一次元的なダムブレーク

これはいけました.ドライベットです.

流速0かつ凸または凹の処理は,中央差分の限界かもしれないですね.

何とか一工夫して使えるようにしたいです.

なお,gifはここにあげるために強烈に画質を落としてます.

f:id:SedimentHydraulics:20161103201324g:plain

水面勾配


集中格子の場合

H(i-1)=H(i+1)=0.0,H(i)=1.0の場合,dH/dx=0 これは、各水面勾配の重みが0.5となり、打ち消しあって0になるようですね。

物理現象からすると、iの地点においてi+1とi-1の地点にそれぞれ流れが発生するようになるはずですが、 集中格子が原因でダメですね。

ちょうど、非圧縮生流体の解析のように、集中格子の場合のチェッカーボードの振動のようですね。

解決策

  • 不等間隔格子

非圧縮生流体の場合、不等間隔格子を使うとチェッカーボードの振動は入らないようですが、 そもそもの解決になっていないように思います。

  • その場しのぎの処理

質量保存がどうなるかわかりませんが、水位をうまく割り振るしかないのかなと思います。 あまりスマートではないですね。

水面勾配項の離散化について

前回の記事で,二次元水面に山のような地形が残る要因がわかりました.

凸地形のdH/dxの処理です.流速が0の場合,中央差分としておりますが,H(i-1)=H(i+1)=0.0,H(i)=1.0の場合,dH/dxは0となってしまいます.このような場合,どういう処理が一般的がご存知ですか?

多分凹地形でもどうようの処理が必要なはずです.