2018-09-30

2次元TMモード電磁波動方程式弱形式その１（とfreefemをはじめた経緯など）

今日は２次元の場合のTMモード電磁波方程式の弱形式について記録しておきます

$\displaystyle \int \Big( { \frac{1}{n^2}\nabla v \cdot \nabla H_z } \Big) d\Omega - \int \Big( { k_0^2vH_z } \Big) d\Omega -\int \Big( { \frac{1}{n^2}\ v \nabla H_z \cdot {\bf n} } \Big) d\Gamma =0$

以下　/* 　から　*/まではコメントなのでお急ぎの方は下へ

/**************************************************************
私がfreefemを取り組み始めたのは
２０１４年の春からだと思う
会社に入って３年目
ちょっとした仕事上の要件で偏微分方程式を解きたい案件に出会ったものの
物理実験系出身で昔、差分法プログラムを実習でちょっと習ったくらいなもの
ぜーんぜん扱ったことありませんでした
そこで積読していた
「有限要素法で学ぶ現象と数理 ―FreeFem++数理思考プログラミング― 」で
上手くいくんじゃね？？と思い出会ったのがきっかけです
その後仕事は大成功でその年、社内で一番大きな技術大会で最年少表彰されました！！

<br />

この本は和書で唯一freefemを解説している本なので一定の価値はあるものの
優れたfreefem独習本であるとはいいがたい出来であると思う

freefemはフリーのソフトなのにすごい！
さらっと実用で使うことさえできてしまう！
みんなももっと使ったらいいし、私ももっと勉強したい！
大学でどんな講義が行われているかはわからないけども（いつもセミナー行きそびれてる。。。）
不慣れな社会人の学習環境が整っているとは言い難い状況だと思います
それの一助になればと＆自分の勉強のためにブログをはじめました

なので多少冗長であっても初心者向けに計算は約さないで書くように努めます
（紙媒体から離れた利点）
ただし見やすさを考慮してまとめは初めに

*****************************************************/

さて改めて
今日は２次元の場合のTMモード電磁波方程式の弱形式について記録しておきます
freefemを扱う場合はまずは自分の扱いたい系の微分方程式（強形式と呼ばれる）から
弱形式と呼ばれる積分型の方程式に直す必要がある

強形式 $\displaystyle \nabla \cdot \Big ( \frac{1}{n^2} \nabla H_z \Big )+ k_0^2 H_z= 0$

座標　ｘ：右正、y：上正、ｚ：紙面手前正
$n$ ：屈折率、 $H_z$ ：z方向の磁場、 $k_0$ ：真空での波数 $(=2\pi/\lambda_0)$ 、 $\lambda_0$ ：真空での波長

これはヘルムホルツ方程式と呼ばれるタイプの方程式で
∇と∇の間に屈折率 $n$ の逆数の二乗が入っているのが印象的な方程式である
屈折率は普通物の境目で、例えばガラスと空気とか ― 不連続で変化するのが物理的に自然な配置
ただ数学的に、そこでの微分∇ってどうなてるんだ？微分方程式として成り立ってんのか？
というようなことを数学徒は考えるようだ

この方程式の導出自体はまたほかの機会に、今日はさらっと弱形式を強形式から導きましょう

強形式 $\displaystyle \nabla \cdot \Big ( \frac{1}{n^2} \nabla H_z \Big ) + k_0^2 H_z= 0$

両辺にテスト関数ｖをかけて積分
(テスト関数はとりあえずはおまじないということにしておきましょうか)
$\displaystyle \int \Big( { v \nabla \cdot \Big ( \frac{1}{n^2} \nabla H_z \Big ) + k_0^2vH_z } \Big) d\Omega= 0$
ここで $d\Omega$ は微小体積素、次元に合わせて体積素、面積素は変えて読む

第一項に第１スカラーのグリーンの定理(部分積分の高次元バージョン)をつかって
First scalar Green's theorem
$\displaystyle \int \Big( { a \nabla \cdot ( u \nabla b) + u (\nabla a)\cdot (\nabla b) } \Big) d\Gamma =\int \Big( { au \nabla b \cdot {\bf n} } \Big) d\Omega$
ここで $d\Gamma$ は微小境界要素、次元に合わせて面積素、辺素は変えて読む

$\displaystyle \int \Big( { \frac{1}{n^2}\ v \nabla H_z \cdot {\bf n}} \Big) d\Gamma -\int \Big( { \frac{1}{n^2}\nabla v \cdot \nabla H_z } \Big) d\Omega$

セットにして

$\displaystyle \int \Big( { \frac{1}{n^2}\nabla v \cdot \nabla H_z } \Big) d\Omega - \int \Big( { k_0^2vH_z } \Big) d\Omega -\int \Big( { \frac{1}{n^2}\ v \nabla H_z \cdot {\bf n} } \Big) d\Gamma =0$

境界条件の積分の部分が重要で
反射壁なのか、吸収壁なのか、波源なのかで
自分で設定していくことになる

書きなれていないのか
ここまで来るのに結構時間がかかったので
今回はここまで

2018-09-23

【注意】複素関数空間の場合の転置コマンド「’」はエルミート転置！

freefem

私みたいにマニュアルに目を通さないで
freefemを扱ってる初心者の方に向けての
注意喚起です

周波数空間でマクスウェル方程式を解きたい場合
転置記号「’」を使いながら内積を使用することになる
例えば、

${ \displaystyle \int{\big( \nabla \times {\bf E}) \cdot (\nabla \times {\bf W}}) d\Omega - k_0^2 \int{ \epsilon {\bf W} \cdot {\bf E} }d\Omega - \int{(n_i \times {\bf W}) \cdot (\nabla \times {\bf E})_i }d\Gamma =0 }$
という数式をコード化すると

Vh<complex> [Ex,Ey,Ez],[Wx,Wy,Wz];

problem FMEQ([Ex,Ey,Ez],[Wx,Wy,Wz])= 
   int3d(Th)(
                Rot(Wx,Wy,Wz)'*Rot(Ex,Ey,Ez)
              - k0^2*([Wx,Wy,Wz]'*[Ex,Ey,Ez])    
             )
 + int2d(Th,In,Out)( 1i*k0*VProduct([Wx,Wy,Wz],[N.x,N.y,N.z])'*VProduct( [Ex,Ey,Ez],N.x,N.y,N.z]) )         
 - int2d(Th,In       )( 2i*k0*VProduct([N.x,N.y,N.z],[Wx,Wy,Wz])'* VProduct([N.x,N.y,N.z],[E0x,E0y,E0z]))
    ;

と書くとする場合がそれにあたる。
このとき「’」はすべてエルミート転置（転置＋複素共役）である。
理解していないと飛び去るエネルギーを計算する際に

complex T = int2d(Th,Out)( 1i*k0*VProduct( [Ex, Ey,Ez],[N.x,N.y,N.z])'*VProduct([Ex,Ey,Ez],[N.x,N.y,N.z]) ) ;

と書くがconj共役をEにつけたくなってしまう
私はこれで半日計算が合わなくて困ったので
皆様は気をつけて

マニュアルの
4.9.3 Matrix Operations
に記述があります
（そのうち本腰入れて
　すみからすみに目を通さないと）

しかしFreeFEMの情報はネット上であんまりないような気がする・・・

2018-09-21

サーモグラフィックのような色諧調

freefem

初めての計算はコードを間違える可能性が大なので
かならずplotで結果のグラフを確認しながらプログラムを行うことが多いです
大体グラフの変になり方からバグの位置が想定されることが多数です

ただfreefemのデフォルトのplotの色諧調はオシャレとはいいがたい。。。
人に見せるときにはそれなりに色諧調をコントロールしたい
結構デザインセンスが要求されるが
よく使う色のセットを作っておくと便利
今回は温度の場合に使いそうな
サーモグラフィックライクな色調を作っておく

サーモグラフの画像からスポイトで色をとってきてhsvに変換
www.peko-step.com

plotコマンドではhsv色空間を線形に補完して
色諧調つないでるっぽい？
ここで問題発生
黄色（ｈ=50/360程度）から赤（ｈ=360/360）に変化させる際に
途中の緑と青をとおって色変化してしまう。。。
そうか！一周できるかもしれない
黄色をh=（360+50)/360で定義してやればOK　（なんだhは1超えても定義できるんだ。。。）

あときれいに繋がったような色身にするのは
結構難しいな
赤のところはもう少し調整必要かも

まあとりあえず
というわけでこんな色身になりました
ほかの色セットもまたつくろ

以下、サンプルコード

//thermographic like color scale & gradation control test 

real[int] xx(10),yy(10);
mesh Th=square(10,10);
fespace Vh(Th,P2);
Vh uh=8*x;

real[int] colorhsv=[ // thermographic color hsv model
410./360.,  0./100. , 100./100., 
410./360., 50./100. , 100./100., 
405./360., 90./100. , 100./100., 
385./360., 85./100. , 100./100.,
365./360., 80./100. ,  85./100.,
345./360., 80./100. ,  65./100.,
325./360., 85./100. ,  50./100.,
300./360., 90./100. ,  35./100.,
280./360., 95./100. ,  20./100.,
260./360.,100./100. ,   0./100.

];

real[int] newviso(81);
for (int i=0;i<newviso.n;i++)
newviso[i]=i*0.1;
plot(uh,value=1,viso=newviso,fill=1,wait=1,hsv=colorhsv);

(20180922追記)
FreeFEMのヴァージョン3.4600000（32bit版）では
ｈが１を超えての定義はできないようで
色が黒くなってしまうので
注意

2018-09-19

freefem並列計算ことはじめ

freefem

これまで結局２次元計算ばかりでかなり良い結果を得てきていたので

なかなか３Dを使ってこなかったのですが

とうとう波動光学も本質的に２Dで近似しにくい系を

扱う羽目になってしまいました

さてとりあえずexampleを動かすことから

http://www.freefem.org/ff++/windows.php

ここによると

MS MPI V7から

msmpisdk.msi

MSMpiSetup.exe　

を落としてどっちもインストール

その後powershellをwindowsのスタートから起動して

cd（チェンジディレクトリ）で

freefemのexample－mpiまで移動して（コマンドライン操作なんてあんまやらない・・・）

とりあえず「mpiexec.exe -np 4 FreeFem++-mpi DDM-Schwarz-Lame-2d.edp」と打ち込む・・・　（うわ、なにこの呪文・・・）

むむ、動いたのか？？

エラーはない

図が出ないとわからん

plotコマンドをコードに入れるもなんも表示されない・・・

コードのコメントを読むと　

「usage :
ff-mpirun [mpi parameter] MPIGMRES2d.edp [-glut ffglut] [-n N] [-k K] [-d D] [-ns] [-gmres [0|1]
argument:
-glut ffglut : to see graphicaly the process」

とあり何やら画面表示したければ「-glut ffglut」を「mpiexec.exe -np 4 FreeFem++-mpi DDM-Schwarz-Lame-2d.edp」の後に続ければよさそう？

そしてじゃじゃん！

おーーー動いた、

詳細はまた解析するとして４つの領域に自動分割してシュワルツ法？で

並列計算させてるのね？

とりあえず最低限完了！！

f:id:Asdfg:20180919040501p:plain

2018-09-10

【エラー報告】3D 周期 missing face

freefem

3次元モデルでフォトニックバンドを計算しようとしている際のエラー
（＊下記はお試しの例）

 missing face 340  146 147 150
 missing face 341  147 150 151
 missing face 342  145 146 149
 missing face 343  146 149 150
 missing face 344  144 145 148
 missing face 345  145 148 149
 3d periodic FE: number of  missing periodic faces 6
  current line = 45
Exec error :  3d periodic condition missing common face
   -- number :1
Exec error :  3d periodic condition missing common face
   -- number :1

どうやら周期関数の関数空間の設定でエラーしているよう
調べてみるとマニュアルNote6.1の特記事項のことらしい

周期的に関数をする場合、周期的にくっつける境界の
四角いメッシュの斜めの線の方向が同じじゃないとエラーをはくが
BuildLayer mesh generatorはメッシュの最大インデックス番号側に
斜めの線を構築するしてしうまうので上手くいかないときがある旨の内容

なので底面のメッシュのナンバリングの再設定が必要
それは以下で行うことができて

 { // for cleanning  memory..
int[int] old2new(0:Th.nv-1);
fespace Vh2(Th,P1);
Vh2 sorder=x+y; 
sort(sorder[],old2new);
int[int]  new2old=old2new^-1;   // inverse the permuation 
//for(int i=0;i< Th.nv;++i) // so by hand. 
//  new2old[old2new[i]]=i;
Th= change(Th,renumv=new2old);
sorder[]=0:Th.nv-1;
}

メッシュに座標の（ｘ＋ｙ）を割り当てて
大きい順に並べなおす感じ？
動作に関しては
examples++-3d/periodic-3d.edpを参照のこと

メッシュの詳細に関して
いじるようなことが初めてだったので
ちょっとわかるまで時間がかかりました
スペシャリストになるには
この辺もわかっておくようにならないとな～～

（BuildLayer mesh generatorのインデックス振りが
　どうなっているのか現状私は理解していないです・・・）

f:id:Asdfg:20180910003247p:plain
f:id:Asdfg:20180910003912p:plain

Asdfgの日記

freefemを使いこなそう

2次元TMモード電磁波動方程式弱形式その１（とfreefemをはじめた経緯など）

*****************************************************/

【注意】複素関数空間の場合の転置コマンド「’」はエルミート転置！

サーモグラフィックのような色諧調

freefem並列計算ことはじめ

【エラー報告】3D 周期 missing face