第60回記念MCMEセミナー

日　時 2023年8月28日（月）13：30～16：40（日本時間）

場　所 武蔵野⼤学有明キャンパス５号館401教室、ハイブリッド開催

FreeFEM - now and in the future (ポスター)

Pierre Jolivet (CNRS) :

“Deep dive into FreeFEM ecosystem”

One of the strengths of FreeFEM is its ability to interact seamlessly with many other scientific libraries, such as MPI for parallel computing, PETSc for (mostly) linear algebra, SLEPc for eigenvalue computation, or HPDDM for domain decomposition methods. In this presentation, I will highlight some design decisions made over the years in order to enable researchers and developers to use FreeFEM as a flexible tool to prototype or implement algorithms, preconditioners,or coupled solvers in different applied fields such as computational fluid dynamics, radiative transfer, solid mechanics.

今回は ICIAM2023 Tokyo で来日された Pierre Jolivet 先生を迎え，大阪大学/理研の鈴木厚先生も加えて，第60回ということで通常より拡大して FreeFEM に関するセミナーを開催いたしました．

FreeFEM は有限要素法を用いた数値計算ツールで，弱形式の方程式をスクリプトで記述することで様々な数理モデルの数値シミュレーションが行えることから，数理工学の分野では広く利用されています．Pierre Jolivet 先生は並列計算ライブラリ HPDDM の開発者であり，その FreeFEM プラグインを実装されています．

Pierre Jolivet 先生からは，FreeFEM におけるスクリプトの基本的な記述について，実際にターミナルから実行して確認しながら説明していただきました．2次元から3次元と切り替えるためのマクロ設定や，ご自身が開発に携わっている PETSc を用いた並列計算の方法などを実演していただきとても参考になりました．

• 講演動画（10分頃から録画）

MCMEseminar-20230828-PJolivet.mp4

※自動作成した英語字幕が見れます（少しづつ誤変換を修正しています）

• Pierre Jolivet さんの web サイト (https://joliv.et/)

  • FreeFEM tutorial  : 今回の講演のスライド

• 高石による introduction

ハンズオン

1. FreeFEM をインストール (できるだけMPI版)

   1. https://freefem.org/

2. サンプルプログラム script.edp を作成

   1. 矩形メッシュの作成

verbosity = 0;

mesh Th = square(20,20);

plot(Th, cmm = "My square");

実行

$ FreeFem++ script.edp

2. 2. 有限要素空間の定義と自由度 (degree of freedom; DOF) の確認

fespace Vh(Th, P2);

cout << Vh.ndof << endl;

Th=square(40,40);

cout << Vh.ndof << endl;

実行

$ FreeFem++ script.edp

2. 3. Poisson  方程式の variational form の定義

Poisson 方程式 : Δu=10,  u=1(矩形領域の下面)

-> 弱形式 (variational form) :  \int_{Th} (du/dx*dv/dx + du/dy*dv/dy + 10*v) dx=0,  u=1 (bottom)

varf vPb(u, v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) + int2d(Th)(10*v) + on(1, u=1.0);

2. 4. 行列と右辺の定義

matrix A = vPb(Vh, Vh);

real[int] rfs = vPb(0, Vh);

2. 5. 変数 u  の定義と求解，表示

Vh u;

u[] = A^-1 * rhs;

plot(u)

実行 (verbosity=0, with-graphics)

$ FreeFem++ script.edp -v 0 -wg

2. 6. 基本境界条件 (essential boundary condition) の計算方法の指定 (tgv=-1 : nonsymmetric elimination)

matrix A = vPb(Vh, Vh, tgv = -1);

real[int] rfs = vPb(0, Vh, tgv = -1);

set(A, solver = GMRES);

実行 (verbosity=0, with-graphics)

$ FreeFem++ script.edp -v 0 -wg

2. 7. IFMACRO を使って，2次元と3次元に対応するプログラムに変更

verbosity = 0;

IFAMCRO(!dimension) // 引数で dimension が与えられない場合（デフォルト）

macro dimenseion() 2 //

ENDIFMACRO

// 2次元

IFMACRO(dimension, 2)

mesh Th = square(20,20);

fespace Vh(Th, P2);

Th=square(40,40);

varf vPb(u, v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) + int2d(Th)(10*v) + on(1, u=1.0);

ENDIFMACRO

// 3次元

IFMACRO(dimension, 3)

load "msh3"

mesh Th = cube(10, 10, 10);

fespace Vh(Th, P1); // 3次元は要素数が多いので P1 メッシュ

Th=cube(5, 5, 5);

varf vPb(u, v) = int3d(Th)(dx(u)*dx(v) + dy(u)*dy(v) + dz(u)*dz(v)) + int3d(Th)(10*v) + on(1, u=1.0);

ENDIFMACRO

// ここから共通部分

plot(Th, cmm="My Mesh");

matrix A = vPb(Vh, Vh, tgv = -1);

real[int] rfs = vPb(0, Vh, tgv = -1);

Vh u;

set(A, solver = GMRES);

u[] = A^-1 * rhs;

plot(u)

$ FreeFem++ script.edp -v 0 -wg

$ FreeFem++ script.edp -v 0 -wg -Ddimension=2

$ FreeFem++ script.edp -v 0 -wg -Ddimension=3

2. 8. メッシュの分割

fespace Vh(Th, P2);

Th = square(40, 40);

plot(Th, cmm = "Before split", wait = 1);

Th = trunc(Th, true, split = 2);

plot(Th, cmm = "After split", wait = 1);

実行 (verbosity=0, with-graphics, dimension=2)

$ FreeFem++ script.edp -v 0 -wg -Ddimension=2

2. 9. submesh の作成と合成

int[int] n2o;

mesh ThSub = trunc(Th, x < 0.5 && y < 0.5, new2old = n2o);

plot(ThSub, cmm = "After submesh extraction", wait = 1);

fespace VhSub(ThSub, P2);

VhSub uSub = cos(y) * sin(4*x);

Vh v;

v = uSub;

plot(ThSub, uSub, cmm = "Finite element function on the extracted mesh", wait = 1);

plot(Th, v, cmm = "Finite element function interpolated on the full mesh", wait = 1);

実行 

$ FreeFem++ script.edp -v 0 -wg -Ddimension=2

2. 10. restriction

int[int] restriction = restrict(Vhsub, Vh, n2o);

v[] = 0;

u = uSub;

v[](restriction) = uSub[]

plot(Th, u, cmm = "Finite element function on the full mesh with connectivity from the parent mesh to the submesh", wait = 1);

実行 

$ FreeFem++ script.edp -v 0 -wg -Ddimension=2

2. 11. 並列計算

cout << mpiSize(mpiCommWorld) << endl; // 各スレッドで並列数を表示

実行 

$ ff-mpi -n 8 script.edp -v 0

2. 12. Overlapping Schwalz methods について

   13. PETSc

load "PETSc"

verbosity = 0;

IFAMCRO(!dimension) // 引数で dimension が与えられない場合（デフォルト）

macro dimenseion() 2 //

ENDIFMACRO

// 2次元

IFMACRO(dimension, 2)

mesh Th = square(20,20);

fespace Vh(Th, P2);

Th=square(40,40);

fespace Vh(Th, P2);

Th = square(40, 40);

plot(Th, cmm = "Before split", wait = 1);

Th = trunc(Th, true, split = 2);

plot(Th, cmm = "After split", wait = 1);

int[int] restriction = restrict(Vhsub, Vh, n2o);

v[] = 0;

u = uSub;

v[](restriction) = uSub[]

plot(Th, u, cmm = "Finite element function on the full mesh with connectivity from the parent mesh to the submesh", wait = 1);

varf vPb(u, v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) + int2d(Th)(10*v) + on(1, u=1.0);

ENDIFMACRO

// 3次元

IFMACRO(dimension, 3)

load "msh3"

mesh Th = cube(10, 10, 10);

fespace Vh(Th, P1); // 3次元は要素数が多いので P1 メッシュ

Th=cube(5, 5, 5);

varf vPb(u, v) = int3d(Th)(dx(u)*dx(v) + dy(u)*dy(v) + dz(u)*dz(v)) + int3d(Th)(10*v) + on(1, u=1.0);

ENDIFMACRO

// ここから共通部分

cout << mpiSize(mpiCommWorld) << endl; // 各スレッドで並列数を表示

plot(Th, cmm="My Mesh");

matrix A = vPb(Vh, Vh, tgv = -1);

real[int] rfs = vPb(0, Vh, tgv = -1);

Vh u;

Mat B(A)

set(B, sparams="-pc_type hypre -ksp_type fgmres -ksp_monitor -ksp_view_final_residual");

u[] = B^-1 * rhs;

plot(u)

実行 (-wg : Graphics 表示)

$ ff-mpi -n 8 script.edp -v 0 -wg

2. 14. 質問