和歌山サイエンスカフェ「インフィニティ」

2017年初頭のモバイルノート（ウルトラブック）な
パソコンを比較検討中。
それなりのCPUパワーと、１６GB、５１２GB（できれば１TB）のSSD
で、１．３キロくらいより軽くて、コンパクトな電源アダプターなのが
必要条件だとする。
お金を出せれば候補も増えるので、２０万円以下くらいが候補。

１．MacBookPro　１３インチ（タッチバーなし）
　　メモリを１６GB、SSDを５１２GBにカスタマイズして、
　　１７９，８００円（税別、学割適用）
　　http://www.apple.com/jp_edu_1460/shop/buy-mac/macbook-pro?product=MLL42J/A&step=config#
　　いまなら、8,500円分のApple Storeギフトカードのプレゼント付き。

　　良い点
　　・英語キーボードが選べる。
　　・Retinaディスプレイになる。
　　・タッチバーなしモデルは、コスパは良好
　　悪い点
　　・CPUがskylake。秋くらいに Kabylakeになるとも言われている。）
・type-C に対応した、SDカードリーダとHDMIソケットのアダプターを買う必要が
　　　自分の場合ある。
　　

２．東芝 dynabook V82/B (東芝ダイレクトだと VZ72/B)
Core i7 7500U で８GBの512GBのSSD(Serial ATA)
２１４、０００円（税別） （東芝ダイレクトの会員になると割引あり）
https://dynabook.com/2in1-mobile-notebook-tablet/v-series/v82b-2016-winter-model-onyx-metallic-12-5-inch-high-end-convertible-pv82bmp-nja.html

　良い点
　・国内メーカーという点
　・いま所有してるR632は、非常に良いPCだった
　・ノングレアな画面で、タッチ対応
　・CPUはKabylake
　・HDMIもVGA接続も可能なアダプター付き（電源供給も可能）

　悪い点
　・英語キーボードなし
　・メモリ８GB設定なし
　・SSDに１TB設定なし。
　・ちょっとベゼル幅が大きめ

３．HP spectre x360（パフォーマンスモデル）
　　CPUはCore i7 7500U、メモリは16GBで、1TBのSSD（NVMe接続）
　　１６１，８００円（税別、学割適用、電話申し込みの場合）
HPの学割サイト　http://jp.ext.hp.com/campaign/personal/notebooks/2017student_pc_collection/
http://jp.ext.hp.com/directplus/personal/?jumpid=ps_74m86jxecn

　良い点
　・CPUはKabylake
　・１６GBも１TBも満たす。しかも、NVMe接続
　・コストパフォーマンつ高め
　・２月２８日までに申し込むと、Office Businessが無料

　悪い点
　・やや重め（１．３ｋｇ）
　・英語キーボードなし
　・電源アダプターのコードが太め
　・HPにやや悪い印象


他にもDELL(XPS13)とかLenovo(Carbon)なども候補にあがるかも。
隠し玉はXiaomiですが、一応、法的な問題もあるので。

spectrex360とXPS13との比較動画


spectre360とMacBookPro13との比較動画

  

第２回目では、さっそくソースをアップする。

今回は、有限体積法の教科書としては、詳しい解説と、例題の丁寧な
回答が掲載されている
『数値流体力学 第2版』
の例題６．２をコード化した。

定常の移流方程式を解いている。
この例題を解いてみたのは、これから解きたい浅水流方程式に
近い形であること。つまり、速度場と密度場が未知であり、
これらを求めるためにSIMPLEアルゴリズムを適用するものであるためである。
また、断面形状も一定でなく、変化しているところも
今後の応用性が高く、参考になると思う。

Javaについては、１か月程度の知識しかなく、
静的メソッドのみの使用と、mainメソッドにほとんどの必要な処理を書いた。
初めて作ったJavaコードなので、至らないところが多い。

勉強しながら、現在はクラス等を用いたコードに修正を加えている。
また、多次元配列の初期化については、
「Arrays#fillメソッドの多次元版の一例」
を利用させていただいた。

連立方程式の解法には、トーマス法を用いている。
格子には等間隔のスタッガード格子を採用している。
また、例題の解答を直接コーディング化ではなく、最後の方に触れられている
不足緩和も係数計算の際に、組み込んでいる。

格子の数や不足緩和の係数を、いろいろ試してみたり、
質量フラックスを計算したりしてみて、有限体積法の面白さを
感じることができると思う。

（注）コードの利用による損害等は一切責任を負いません。

/*
* 有限体積法の計算
* 数値流体力学（第２版）
*
* 2017/2/11
* KZ
*
*/

package FVM;

import java.util.Arrays;

public class FvmTest007 {

public static void main(String args[]){

int nn = 5; // 圧力ポイントの数
int it_max = 200; // 最大回数
int count = 0; // 計算回数の記憶用
double r_limit = 10e-5; // 残差の制限


double m_0 = 1.0; //　初期仮定の質量フラックス
double p_0 = 10.0; // 上流端の圧力
double rho = 1.0; // 密度
double alpha = 0.7; // 不足緩和の係数 (nn,alpha)=(10,0.7),(20,0.5)
double residual = 0.0; // 残差

//double[] A_p = {0.5 ,0.4 ,0.3 ,0.2 ,0.1}; //圧力ポイントの面積
double[] A_p = new double[nn];
for(int i=0; i < nn; i++){
A_p[i] = 0.5 - (0.5-0.1)/(nn-1) * i;
}
//double[] A_u = {0.45,0.35,0.25,0.15};
double[] A_u = new double[nn-1];
for(int i=0; i < nn-1; i++){
A_u[i] = ( A_p[i] + A_p[i+1] ) / 2.0;
}

double[] x_u = new double[nn-1]; // 速度場
double[] u_star = new double[nn-1]; // 補正後の速度
double[] u_old = new double[nn-1]; // 速度の古いもの保存（残差計算用）
//double[] x_p = new double[nn]; // 圧力場
//double[] x_p = {10.0, 7.5, 5.0, 2.5, 0.0}; // 圧力場と仮定初期値
double[] x_p = new double[nn];
for(int i = 0; i < nn; i++){
x_p[i] = 10.0 - (10.0 - 0.0)/(nn-1) * i;
}
double[] p_prime = new double[nn]; //　圧力の補正値
double[] p_star = new double[nn]; // 補正後の圧力
double[] d = new double[nn-1]; // 圧力補正のパラメータ

double u_in = 2.0; // 上流端の流速（境界条件）

//方程式の係数
double F_w = 0.0;
double F_e = 0.0;
double a_W = 0.0;
double a_E = 0.0;
double a_P = 0.0;
double S_u = 0.0;

double b_prime = 0.0;

//連立方程式の係数行列
double[][] matrixA = new double[nn-1][nn-1];
double[] matrixB = new double[nn-1];


double[][] matrixC = new double[nn][nn];
double[] matrixD = new double[nn];


//配列の初期値
Arrays.fill(x_u,0.0);
Arrays.fill(u_star, 0.0);
Arrays.fill(p_prime, 0.0);
Arrays.fill(p_star, 0.0);
Arrays.fill(d, 0.0);
ArraysfillEx.fill(matrixA, 0.0); //多次元配列の初期値
Arrays.fill(matrixB, 0.0);
ArraysfillEx.fill(matrixC, 0.0); //多次元配列の初期値
Arrays.fill(matrixD, 0.0);



//速度場の初期仮定の設定
for(int i = 0; i < nn-1 ; i++){
x_u[i] = m_0 / ( rho * A_u[i] );
}

System.out.println("x_u : "+Arrays.toString(x_u)); // 行列の標準出力
System.out.println("x_p : "+Arrays.toString(x_p)); // 行列の標準出力
System.out.println("A_u : "+Arrays.toString(A_u));
System.out.println("A_p : "+Arrays.toString(A_p));



//計算の開始
for(int it = 0; it < it_max ; it++){

System.out.println("Count : " + it);

// 速度の更新（残差計算用）
for(int i = 0; i < nn-1 ; i++){
u_old[i] = x_u[i];
}
residual = 0.0; // 残差の初期化


//運動量保存則　始まり
//速度の格子点//
//格子点　1(0)
u_in = x_u[0] * (A_u[0] / A_p[0]); System.out.println("u_in : " + u_in);
F_w = rho * u_in * A_p[0]; //System.out.println(F_w);
F_e = rho * (x_u[0] + x_u[1])/2 * A_p[1]; //System.out.println(F_e);

a_W = 0.0;
a_E = 0.0;
a_P = F_e + F_w * 0.5 * Math.pow((A_u[0]/A_p[0]),2.);


S_u = (x_p[0]-x_p[1])*A_u[0] + F_w*(A_u[0]/A_p[0])*x_u[0] + (1-alpha)*a_P/alpha * x_u[0];
a_P = a_P/alpha; //不足緩和

matrixA[0][0] = a_P;
matrixA[0][1] = -1.0 * a_E;

matrixB[0] = S_u;

d[0] = A_u[0]/ a_P * alpha;//不足緩和


//格子点　2,3(1 〜 nn-3)
for(int i=1 ; i < nn-2;i++){
F_w = rho * (x_u[i-1] + x_u[i])/2 * A_p[i];
F_e = rho * (x_u[i] + x_u[i+1])/2 * A_p[i+1];
a_W = Max(F_w,0.0);
a_E = Max(0.0, -1.0*F_e);
a_P = a_W + a_E + (F_e - F_w);

S_u = (x_p[i]-x_p[i+1])*A_u[i] + (1-alpha)*a_P/alpha * x_u[i];//不足緩和
a_P = a_P/alpha; //不足緩和

matrixA[i][i-1] = -1.0 * a_W;
matrixA[i][i] = a_P;
matrixA[i][i+1] = -1.0 * a_E;

matrixB[i] = S_u;

d[i] = A_u[i]/ a_P * alpha;
}


//格子点 4(nn-2)
F_w = rho * (x_u[nn-3] + x_u[nn-2])/2 * A_p[nn-2]; //System.out.println(F_w);
//F_e = rho * (x_u[nn-2] + x_u[1])/2 * A_p[1]; //System.out.println(F_e);
if(it == 0) {
F_e = m_0; // 最初の計算だけ　要検討！！！
}else{
F_e = rho * x_u[nn-2] * A_u[nn-2];
}

a_W = F_w;
a_E = 0.0;
a_P = a_W + a_E + (F_e - F_w);

S_u = (x_p[nn-2]- x_p[nn-1])*A_u[nn-2] + (1-alpha)*a_P/alpha * x_u[nn-2];
a_P = a_P/alpha; //不足緩和

matrixA[nn-2][nn-3] = -1.0 * a_W;
matrixA[nn-2][nn-2] = a_P;

matrixB[nn-2] = S_u;

d[nn-2] = A_u[nn-2]/ a_P * alpha;



//System.out.println("matA: "+Arrays.deepToString(matrixA));
//System.out.println("matB: "+Arrays.toString(matrixB));
//System.out.println(" d : "+Arrays.toString(d));

//連立方程式を解く

u_star = thomas(nn-1,matrixA,matrixB);

System.out.println("u_star: " + Arrays.toString(u_star));

//運動量保存則　終わり

//圧力の補正式　始まり
//圧力の格子点 1,A (0)
matrixC[0][0] = -1.0;

//圧力の格子点 2から4,A (1 - nn-2)
for(int i = 1; i < nn-1 ; i++){
a_W = rho * d[i-1] * A_u[i-1];
a_E = rho * d[i] * A_u[i];
F_w = rho * u_star[i-1] * A_u[i-1];
F_e = rho * u_star[i] * A_u[i];
a_P = a_W + a_E;
b_prime = F_w - F_e;

matrixC[i][i-1] = -1.0 * a_W;
matrixC[i][i] = a_P;
matrixC[i][i+1] = -1.0 * a_E;

matrixD[i] = b_prime;
}

//圧力の格子点 5,E (nn-1)
matrixC[nn-1][nn-1] = -1.0;


//System.out.println("matC: "+Arrays.deepToString(matrixC));
//System.out.println("matD: "+Arrays.toString(matrixD));

//連立方程式を解く
p_prime = thomas(nn,matrixC,matrixD);

System.out.println("p_prime: " + Arrays.toString(p_prime));

//圧力の補正式　終わり

//圧力の補正　始り
for(int i=0;i < nn;i++){
x_p[i] = x_p[i] + alpha * p_prime[i];
}
System.out.println("calculated x_p : " + Arrays.toString(x_p));
//圧力の補正　終わり

//速度の補正　始まり
for(int i=0;i < nn-1;i++){
x_u[i] = u_star[i] + d[i] * (p_prime[i] - p_prime[i+1]);
}
System.out.println("calculated x_u : " + Arrays.toString(x_u));
//速度の補正　終わり
//格子点Aの圧力の補正
x_p[0] = p_0 - 0.5 * rho * Math.pow(x_u[0], 2.0) * Math.pow(A_u[0]/A_p[0],2.0);
System.out.println("calculated2 x_p : " + Arrays.toString(x_p));
//　終わり

// 残差の計算
for(int i = 0; i < nn-1 ; i++){
residual += Math.pow(u_old[i]-x_u[i],2.0);
}
System.out.println("residual : " + residual);

if(residual < r_limit) {
count = it;
break;
}

}

System.out.println("¥n Count: " + count);
System.out.println("x_u : " + Arrays.toString(x_u));
System.out.println("x_p : " + Arrays.toString(x_p));
System.out.println("A_u : "+Arrays.toString(A_u));
System.out.println("A_p : "+Arrays.toString(A_p));



}



//以下、メソッド

//大きい数字を返す
public static double Max(double a,double b){
if(a>b){
return a;
} else {
return b;
}
}


//トーマス法
public static double[] thomas(int n,double[][] A,double[] B){

double[] x = new double[n];

double[] g = new double[n];
double[] s = new double[n];

//配列の初期値
Arrays.fill(x,0.0);
Arrays.fill(g, 0.0);
Arrays.fill(s, 0.0);

g[0] = A[0][0];
s[0] = B[0];

for(int i = 1 ; i < n ; i++){
s[i] = B[i] - A[i][i-1]*s[i-1]/g[i-1];
g[i] = A[i][i] - A[i][i-1]*A[i-1][i]/g[i-1];
}

x[n-1] = s[n-1]/g[n-1];

for(int i = n-2 ; i >= 0; i--){
x[i] = (s[i]-A[i][i+1]*x[i+1])/g[i];
}

return x;
}
}

  

windows7ノートに
Latexをインストールした。

Texインストーラ３を利用
http://www.math.sci.hokudai.ac.jp/~abenori/soft/abtexinst.html

参考にしたサイト
http://did2memo.net/2016/04/24/easy-latex-install-windows-10-2016-04/  

ここ数ヶ月、流体の数値計算の勉強を続けている。
２次元の河床変動プログラムを作るのが、当面の最終目標。

河床変動の計算の基礎式は、
１：流体の質量保存則（連続の式）
２；流体の運動量保存則（ナビエストークス方程式）
３；流砂の連続式
４；流砂量式
である。これら４つの方程式を解くことで
流れの場を一般的には求めることができる。

これらの式は複雑に絡み合い、特に２の運動量保存則は、
非線形の方程式（移流方程式）であり、一般に解析的に解くことが
できない。
そこで、数値計算が多用される。

数値計算においては、上記の方程式を離散化して、
解く方法がとられている。

その離散化の方法の違いにより、大きく
１：　格子法
２：　粒子法
に区別される
（他に個別要素法があるが、どう分類されるのかはまだ、勉強不足）

最近は２の粒子法も注目を集めているが、
計算範囲が大きい場合や、時間スケールを長くするには、
計算パワーが必要なことから、まだまだ一般的ではない。
粒子法については
『粒子法入門 流体シミュレーションの基礎から並列計算と可視化まで C/C++ソースコード付』
が詳しくて分かりやすい。

このため、一般には１の格子法が採用される。
格子法は、さらに、
A：差分法
B：有限体積法
C；有限要素法
などに分類される。

差分法は、微分方程式をテイラー展開を利用して離散化する。
また、有限要素法は、重み付きの関数を用いて離散化する。
（有限要素法も勉強不足）
一方で、有限体積法は、コントロールボリュームと呼ばれる
有限な体積の要素に積分することで、その収支を用いて離散化する。
有限体積法では、コントロールボリュームの質量収支等を考慮するので、
質量が保存されやすいという特があり、この点が他の計算方法との違いと
言われている。

連続体としての質量の保存を重視する流体力学の数値計算では、
有限体積法が人気を集めているようである。
他の工学分野（構造力学等）では、有限要素法が一般的なようである。

以上のような理由から、有限体積法を用いて、
数値計算（数値シミュレーション）を進めていくこととした。


なお、使用言語は、
Javaを使っていくことにする（途中で変更するかもしれない）
もともと数値計算でよく用いられているFortranも使っていたが、
Androidアプリ開発のため、Javaを勉強し始めたことと、
オブジェクト指向を理解するためにも、Javaをプログラミング言語として
学んでいくことは、今後に役立つと考えたためである。


有限体積法の勉強に用いた参考文献を以下に示す。