気になるなあ。

プロジェクトページに論文もありました。

http://www.byungsoo.me/project/deep-fluids/paper.pdf

論文とオープンソースのコラボは最強だな。憧れますね。

2020-01-27

浮遊砂の濃度分布：Rouse分布と板倉・岸の式

久しぶりに浮遊砂のプログラムを書いているのでチェック用に書いてみました。

Rouse分布で横軸を濃度に変更
基準面濃度式は板倉・岸の式
沈降速度はRubeyの式
河床勾配1/1000、水深5mのケース

ソースは以下です。

Jupyter Notebook Viewer

たかだかこれだけなのに式は超煩雑。
こう見るとWashloadって0.1mmくらいなんでしょうね。
そう言えば、土砂水理学1の板倉・岸の式（5.3.95）が間違ってました。

リンク

2020-01-24

備忘録：pdfをsvgに変換する

いつもInkscape使ってたけど、WSL-ubuntuだとこれで一発。複数ページは無理だけど。

 pdftocairo -svg test.pdf test.svg

Ubuntu Manpage: pdftocairo - Portable Document Format (PDF) to PNG/JPEG/TIFF/PDF/PS/EPS/SVG using cairo

シェルスクリプトでの連続実行はこんな感じ

#!/bin/sh
files="*.pdf"
for filepath in ${files}
do
  echo $filepath
  pdftocairo -svg "${filepath}" "${filepath}.svg"
done

なお、インストール方法は以下のとおりです。

sudo apt update

sudo apt -y install poppler-utils

初めてgoogle Colaboratoryを使ってみましたがやばいですね。いきなりtensorflowが使える。ちょっとした作業ならanacondaより圧倒的に楽。

colab.research.google.com

2020-01-17

多列砂州河道の再現計算（失敗例）

はじめに

こんな感じの多列砂州形成の計算をしてみたい。

text

静岡河川事務所HPよりリンク

計算結果（失敗例）

計算条件

平面二次元流+掃流砂単一粒径河床変動
流下方向13km × 川幅650m
格子サイズは25m x 25m
粒径40mm
周期境界条件
初期条件：流れ水深1.6m 河床：平均河床勾配1/400に対して粒径サイズの凹凸をランダムに与える。
計算時間：20000時間(25日間にかかりました）

計算結果

8000時間後までの河床変動量と水深を図化
最初のほうは多列砂州が形成されるが、時間の経過とともに単列化する。最終的にはほどんど動かない。
大井川等では過去100年以上は多列砂州が維持されている。つまり、計算がおかしい。

f:id:SedimentHydraulics:20200117171644g:plain

考察

計算では河床の擾乱が消えていくが、実現象ではそうではないため、擾乱が生成される要因が存在する。
黒木・岸の研究（https://www.jstage.jst.go.jp/article/jscej1969/1984/342/1984_342_87/_article/-char/ja/）では、その要因は流砂の遅れ時間（非平衡性）と述べられている。
近年若い先生方には、実験・解析より多列砂州は維持されないという研究成果がいくつか見られるがそれは説明がつかないと思う。
私は局所勾配の流砂量への影響が効くのではと考えており、それを考慮した計算をやってみようと考えている。

その他

流体系？の先生の論文（http://traffic.phys.cs.is.nagoya-u.ac.jp/~mstf/pdf/mstf2003-01.pdf）をたまたま読んだのですが、浮遊土砂の拡散がないと河岸ができないと書いていました。河川系では、流路形成の基本は掃流砂ということは常識なので違和感が半端ないです。実際に私の計算でも掃流砂のみ流路、河岸はできているし。

2020-01-12

numba.jit解説1 : 単純な反復計算の高速化

このブロクの人気記事であるnumba.jit関連について、もう少し詳しく書いてみました。

computational-sediment-hyd.hatenablog.jp

私はnumba.jitのヘビーユーザーです。おかげで数日かかるような科学技術計算もpythonで書くようになりました。

numba.jitでは単純な処理の反復が高速化されます。そのため、幾何演算や科学技術計算で効果的です。

自分の勉強のためにも、numba.jitの特徴をまとめておきます。

環境

OS:windows10 pro, CPU Core i7-8550U
jupyter上で動作確認、時間計測は%%timeitを使用
numba 0.47

簡単な事例で動作確認：ポリゴンの面積の計算

ソースコード

ポリゴンの面積を計算するコードでnumba.jitを解説します。

面積を計算する式はこんな感じです。

$S = \frac{1}{2} \left| \sum^n_{i=2}{(x_{i-1}y_{i} - x_{i}y_{i-1})} \right|$

shapelyで半径1の円を簡素化して多角形を作成。64角形になりました。

from shapely.geometry import Point

p = Point(0.0, 0.0)
x = p.buffer(1.0)
s = x.simplify(0.001, preserve_topology=False)
polygon = np.array( [ p for p in s.exterior.coords])
len(polygon)
# 65

f:id:SedimentHydraulics:20200112134907p:plain

上のポリゴンの面積を求めるプログラムはこんな感じでしょうか。当然面積はほぼ3.14になります。

import numpy as np

def poly(p):
    s = np.array( [ p[i-1][0]*p[i][1] - p[i][0]*p[i-1][1] for i in range(1, len(p)) ] )
    return 0.5*np.abs(np.sum(s))

poly(np.array(polygon))

numba.jitを使う場合は、defの前に@jit(nopython=True)を付けて、

import numpy as np
from numba import jit

@jit(nopython=True)
def poly_numba(p):
    s = np.array( [ p[i-1][0]*p[i][1] - p[i][0]*p[i-1][1] for i in range(1, len(p)) ] )
    return 0.5*np.abs(np.sum(s))

です。ほとんど何も変わりません。

速度計測

さっそく上の2つのコードの速度比較です。

%%timeit
poly(np.array(polygon))

105 µs ± 1.04 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%%timeit
poly_numba(np.array(polygon))

1.1 µs ± 34.8 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

なんとこれだけで約100倍の速度差が出ます。

注意点:いつでも使えば良いわけではない。

numba.jitは関数の定義時に通常の場合と比較して若干時間がかかります。関数の定義を含めて時間計測を行うと、

%%timeit 

def poly(p):
    s = np.array( [ p[i-1][0]*p[i][1] - p[i][0]*p[i-1][1] for i in range(1, len(p)) ] )
    return 0.5*np.abs(np.sum(s))

poly(np.array(polygon))

105 µs ± 931 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%%timeit

@jit(nopython=True)
def poly_numba(p):
    s = np.array( [ p[i-1][0]*p[i][1] - p[i][0]*p[i-1][1] for i in range(1, len(p)) ] )
    return 0.5*np.abs(np.sum(s))

poly_numba(np.array(polygon))

151 ms ± 3.31 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

のとおり、通常の関数の方が早くなります。そのため、一回だけ関数を呼び出す場合は不利になります。この例では10万回呼び出してようやくnumba.jitのほうが早くなります。

デメリット：pythonらしい書き方が使えない。

今回の例をpythonらしく書くと例えばzip関数を使って、

def poly(p):
    s = np.array( [ p0[0]*p1[1] - p1[0]*p0[1] for p1, p0 in zip(p[1:], p[:-1]) ] ) 
    return 0.5*np.abs(np.sum(s))

と書けますがnumba.jitではエラーがでます。

numbaのupdateで徐々に使える関数が増えてますが他にもnumpyやscipyの関数も使えないものが多いです。

まとめ

numba.jitは大規模演算用
大規模演算では必須
次回はクラス編を書きます。

computational-sediment-hyd.hatenablog.jp

2020-01-09

epsg.ioが壊れている？

2020/2/1追記

いつの間にか修正されていました。もう普通に使えます。

Endpoints now throwing 500 errors · Issue #156 · maptiler/epsg.io · GitHub

http://epsg.io/

のAPIが昨日から壊れているみたいです。

githubにはいくつかissuesが上がってました。

https://github.com/SciTools/cartopy/issues/1437

https://github.com/rhattersley/pyepsg/issues/15

matplolib,holoviews,geoviews等pythonの地図系のモジュールの多くで使われるcartopyモジュール内ではpyepsgモジュールを参照していますが、この中でepsg.ioのAPIが使われてます。

そのため、pythonでepsgを使った地図のレンダリングが全く使えません。

結構大事かなと思ったのですがほとんど情報がないです。

なんとかならないかな。。。

2020-01-06

本年もよろしくお願いいたします

本年もよろしくお願いいたします。

本業が忙しくて全然書けおりませんが今年は更新頻度を上げたいです。

このブログ（勉強会）での今年中の目標は以下の3点。

小規模河床波上の流れのモデルの開発：鉛直二次元流
網状流路河道の自律的河床変動機構モデルの開発：平面二次元流
地域防災に資するwebアプリ開発

引き続きよろしくお願いいたします。

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

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

Deep Fluidsが気になる