趣味で計算流砂水理

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

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

MENU

お久しぶりです

いろいろ書けなくてすみません。

最適化数学

最適化数学の場合、以下の書籍を読んでました(過去形)。

入門者の本です。

http://www.kyoritsu-pub.co.jp/bookdetail/9784320017863

ただどちらかというと、線形計画のみピックアップした読み方をしてたので、 無制約の部分が使えるかどうかは不明です。

計算時間がクリティカルって面白いですね。 Fortranでやってもかなり計算がかかるものなのでしょうか?

Fortranで問題なければ、エクセルからDLLをキックすれば行けそうな気がするのですが。

グラフ

  • すごいです。綺麗に収束してますね(感動)。
  • 中のアルゴリズムが気になりますが、面白い研究で実用も使えて良いものだと思います。
  • なんかちゃんと理論立てて、数学を使えばちゃんと上手くやれるんですね。当たり前ですが・・・。

Flux Upwind Scheme

  • 2次元版を書いてみましたが、振動が入ってるっぽくて、どこか間違ってるようです。
  • 今度教えてください。
  • Parabolic Bowl でも試してみます。

f:id:SedimentHydraulics:20160911110429p:plain

VR

  • VRよりはARが流行りそうかなと思ってます。
  • ここら辺だと、ゲーム関連をやってる人だと詳しそうです。

Deep Learning

  • どうも、深層学習は画像しか上手くいっていないみたい、な話が出ました。
  • 画像は、JPGであれば、フーリエ変換を応用したものを使っており、親和性がもともと高いのではないかということです。

無制約最適化問題

ぜんぜん書けなくてすいません。完全に執筆モードです。

無制約最適化問題

ここのところ、夢中になっているアジョイント法ですが、ポイントの一つに降下法の解法があります。 基本は、共役勾配法(正確には、非線形共役勾配法)を使ってますが、ステップ幅の問題があり、準ニュートン法(正確には、記憶制限準ニュートン法)にしたいなと思ってます。

このような話って、無制約最適化問題の反復法による解法に当たると思いますが、良い教科書的な本ってご存知ですか? なかなか見つからなくて困ってます。

ジョイント法ですが、実河川スケールに拡張した場合、行列が大きくなり過ぎて、降下法の演算に時間がかかり過ぎます。なんとかしないとまずい感じです。

勉強会のテーマについて

ちょっと本業がバタバタしており、あまり準備ができそうにないのですが、私からは以下のような準備をしています。

  • 平面二次元の話

    • Upwind Flux Schemesの境界条件を教えて下さい。
    • Upwind Flux Schemesの離散化について
    • parabolic bowlの話
  • 久しぶりに断面準二次元の話

    • 今後の展開について私から報告です。
  • その他

こんな感じでしょうか。他に思いついたら準備しておきます。

Linuxとか

書き込み遅くなりました。

Linux

サーバーで動かすという場合であれば、CentOSが多いような気がします。 WebもCentOS(もしくはFree BSD)が多いようです。

チェック用であれば、VMwareでも良いかと思いますが、 実際に稼働するサーバーと同じ環境にしておいたほうが良いかと思います。

また、ライブラリを利用しているようでしたら、そこらへんにも気を使う必要があるかと思います。

個人的にはCentOSをネイティブに用意しておくのも悪くないと思います。

なお、個人PCにCent OS 6は入れてます(Win7とのデュアルブート)が、最近使ってません・・・。

Parabolic bowl

これでチェックできそうですね。 まだ回せていませんが、Python良いです。 リスト内包表記は強烈です(笑)

中身は後ほど確認します。

Pythonでプロトタイプを作って回しますが、 結局中身が肥大化して遅くなるので、結局Fortranに書き換えたくなります。。。 この前教えてもらったFortranでのDLLを試してみたいと思います。

最近

いろいろあって、波動方程式とかシェルスクリプトとかPythonとか仕事で使いそうです。 まあ、WinなのでCygwinですが・・・。

Parabolic Bowl

更新空いてしまいました。

Parabolic Bowl

先週末から平面二次元の計算を始めました。 手始めにparabolic bowlの図化をpythonで書いてみました。 気が向いたら回してみてください。DBに置いておきます。 さすがに遅すぎるので実計算には使えそうにないです。動画っぽくしたいなら、ImageMagickでgifにするほうが良いかもしれないです。

f:id:SedimentHydraulics:20160713185217p:plain

ソースコードはこんな感じです。ほぼグラフのところで、解析解の計算は数行です。やっぱりI/Oは手間ですね。

# -*- coding: UTF-8 -*-
# use python3.4
## @package
#
#
#
###%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np

class MyGragh:
    def __init__(self,X,Y,zb,wL,idX,idY):
        self.X,self.Y,self.idX,self.idY = X,Y,idX,idY
    #make gragh
        self.fig = plt.figure(figsize=(12, 10))
    #make 3Dgragh
        self.ax = self.fig.add_subplot(2, 1, 1, projection='3d')
        self.ax.set_zlim3d(0, 14)
        self.ax.view_init(elev=33, azim=38)
        Xm, Ym = np.meshgrid(self.X, self.Y)
        h = list(map(list,zip(*zb)))
        self.ax.plot_surface(Xm, Ym, h, color='k', zorder=1, rstride=1, cstride=1,linewidth=0.1, antialiased=False,alpha=0.3)
        h = list(map(list,zip(*wL)))
        self.surf = self.ax.plot_surface(Xm, Ym, h, cmap=cm.jet, zorder=2, rstride=1, cstride=1,linewidth=0.1, antialiased=False, alpha=1)
    #make 2Dgragh
        self.ax2 = self.fig.add_subplot(2, 2, 3)
        plt.title('X = '+str(self.X[self.idX]))
        h = zb[self.idX,:]
        plt.plot(self.Y,h, color='k')
        h = wL[self.idX,:]
        self.lineX = self.ax2.plot(self.Y,h, color='b')

        self.ax3 = self.fig.add_subplot(2, 2, 4)
        plt.title('Y = '+str(self.Y[self.idY]))
        h = zb[:,self.idY]
        plt.plot(self.X,h, color='k')
        h = wL[:,self.idX]
        self.lineY = self.ax3.plot(self.X,h, color='b')

    def plotWL(self,wL,t,u,v):
    #make 3Dgragh
        self.surf.remove()
        Xm, Ym = np.meshgrid(self.X, self.Y)
        h = list(map(list,zip(*wL)))
        self.surf = self.ax.plot_surface(Xm, Ym, h, cmap=cm.jet, zorder=2, rstride=1, cstride=1,linewidth=0.1, antialiased=False, alpha=1)
    #make 2Dgragh
        for c in self.lineX : c.remove()
        h = wL[self.idX,]
        self.lineX = self.ax2.plot(self.Y,h, color='b')
        for c in self.lineY : c.remove()
        h = wL[:,self.idY]
        self.lineY = self.ax3.plot(self.X,h, color='b')

    #make title
        self.fig.suptitle('t='+str(t)+'[sec]'+' U='+str(round(u,3))+'[m/s]'+' V='+str(round(v,3))+'[m/s]')
        plt.pause(0.001)

    def __enter__(self):
        return self

    def __exit__(self, exc_type, exc_value, traceback):
        plt.savefig('output.png')
        return False

# main
nx, ny = 21, 21
X = np.arange(-10, 11, 1)
Y = np.arange(-10, 11, 1)
a = 7
h0 = 4
zb = np.array([[h0/a**2*(X[i]**2+Y[j]**2) for j in range(ny)] for i in range(nx)])
eta = a/3
tPriod=60
omega = 2*np.pi/tPriod
#wL = lambda t: np.array([[ eta*h0/a**2*(2*X[i]*np.cos(t*omega)+2*Y[j]*np.sin(t*omega)-eta)+h0 for j in range(ny)] for i in range(nx)])
wL = lambda t, x, y: eta*h0/a**2*(2*x*np.cos(t*omega)+2*y*np.sin(t*omega)-eta)+h0
wLr = lambda t: np.array([[ wL(t,X[i],Y[j]) if wL(t,X[i],Y[j]) - zb[i,j] > 0 else zb[i,j] for j in range(ny)] for i in range(nx)])
u = lambda t: -eta*omega*np.sin(t*omega)
v = lambda t: -eta*omega*np.cos(t*omega)

idX=10
idY=10
tmax = 600
with MyGragh(X,Y,zb,wLr(0),idX,idY) as g:
    for t in range(1,tmax):
        if t%10==0 : g.plotWL(wLr(t),t,u(t),v(t))

ジョイント

とりあえず、矩形断面ではできたっぽくて、テンションが上がってます。パッと見せれる絵がないので勉強会のときにでも話します。

linux環境を構築したい

linux環境を構築したいわけ

ずっとlinux環境が良いと考えておりましたが、以下の理由によりそろそろやろうかと

OSに依存しないコードの作成

このブログで書いているコードをある程度できたものはgithubで管理しようかなと思ってます。でも、例えば、先日のpython + fortran DLLのコードとかはwindowsに限定しているような気がします。(configureとかをしっかり書けば良いのですが。。。)そんなこともあり、開発環境としてwindowsはどうかなとようやく思い始めました。

諸事情でCentOSで動くコードを作る必要がある

仕事でちょっとした水理計算をクラウド上のCentOSで動かす必要があり、チェック用にlinux環境が欲しいなと。

質問

VMware or デュアルブート ...etc

  • できればモバイルで作りたいです。そうなるとwindowsベースでVMware or デュアルブートかなと思ってます。
  • それとも新たなモバイルを買って作ったほうが良いでしょうか。

おすすめのlinuxディストリビューション

  • 基本的なところですとCentOSとかUbuntuとかでしょうか。いろいろあって面白そうです。思い切って軽量なelementary OSとかでもいいですが。
  • せっかくなのでいっしょに入れませんか。

その他

  • 2D SWEの話はいろいろ書きたいので後日。
  • 最近EUROの見過ぎでいろんなことが遅れ気味です。