分散離散化でマルチグリッド前処理行列を使用して微分方程式を解く

ライブスクリプトを開く

この例では、前提条件を使った反復ソルバーと分散配列を使用してポアソン方程式を解く方法を示しています。分散配列を使用すると、1 台のマシンのメモリだけでなく、マシンのクラスターのメモリも使用して、計算をスケールアップできます。スケールアップはコードを変更せずに実行できます。

この例は、分散配列を使用した反復法による線形方程式系の求解で説明されているトピックの続きです。[1] に基づいて、この例はポアソン方程式を使用して室内の熱分散を、均質定常熱伝導方程式として知られる形式でモデル化します。定常とは時間が経過しても熱は変化しないことを表し、均質とは外部熱源がないことを表します。

$- Δ u = - (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} + \frac{\partial^{2} u}{\partial z^{2}}) = 0$

この方程式の $u$ は室内の各点 $(x, y, z)$ における温度を表します。この方程式を解くには、まず有限差分離散化法を使用して、線形方程式系により近似します。その後、前処理付き共役勾配法 (pcg) を使用して、この方程式系を解きます。前提条件により問題が変換され、数値ソルバーのパフォーマンスが向上します。分散配列を使用すると、マシンのクラスターの結合メモリを利用して、より精密な離散化を行うことができます。

以下の方法について説明します。

離散 3 次元グリッドと境界条件を設定する。
離散化とマルチグリッド前処理行列を定義する。
前処理を使用した数値ソルバーを適用して、3 次元ボリュームにまたがる熱伝導方程式を解く。

空間次元の離散化

この例では、一片が 1 の立方体で部屋をモデル化しています。最初に、3 次元グリッドを使用してこれを離散化します。

この例での前処理法では、粒度のレベルが異なる複数のグリッドを使用しています。各レベルで、次元ごとに係数 2 を使用してグリッドが疎化されます。マルチグリッドレベルの数を定義します。

multigridLevels = 2;

最も粒度が高いグリッドの X、Y、Z の各次元の点の数を定義します。この前処理法では、このグリッドの点の数が 2^multigridLevels で割り切れることが必要です。この例では、マルチグリッドレベルの数が 2 なので、点の数は 4 で割り切れなければなりません。

numPoints.X = 32;
numPoints.Y = 32;
numPoints.Z = 32;

関数 meshgrid を使用して、3 次元グリッドにより空間次元を離散化します。linspace を使用して、点の数に従って各次元を均等に分割します。立方体の境界が含まれるようにするには、2 つの点を追加しなければなりません。

[X,Y,Z] = meshgrid(linspace(0,1,numPoints.X+2), ...
    linspace(0,1,numPoints.Y+2), ...
    linspace(0,1,numPoints.Z+2));

境界条件の定義

部屋に窓とドアがあるものとします。壁と天井の温度は 0 度で一定で、窓の温度は 16 度で一定で、ドアの温度は 15 度で一定です。床は 0.5 度です。目標は、室内全体での温度分布を判別することです。

床、窓、ドアの座標を関係演算子を使用して定義し、これらの境界要素の温度を定義します。境界は立方体の小平面です。そのため、X、Y、Z のいずれかが 0 または 1 でなければなりません。残りの境界と立方体の内部は 0 に設定します。

floor  = (0.0 <= X & X <= 1.0) & (0.0 <= Y & Y <= 1)   & (Z == 0.0);
window = (X == 1)              & (0.2 <= Y & Y <= 0.8) & (0.4 <= Z & Z <= 0.6);
door   = (0.4 <= X & X <= 0.6) & (Y == 1.0)            & (0.0 <= Z & Z <= 0.6);

u = zeros(size(X));
u(floor)  = 0.5;
u(window) = 16;
u(door)   = 15;

これらの境界条件では、領域の境界に沿って解がとらなければならない定数値が指定されます。このような種類の境界条件はディリクレ境界条件と呼ばれています。

関数 slice を使用して境界条件を可視化します。立方体の境界に配置されたスライスを使用し、非ゼロの境界条件を表示します。

xSlices = 1;
ySlices = 1;
zSlices = 0;
f = slice(X,Y,Z,u,xSlices,ySlices,zSlices,'nearest');
title('Constant nonzero boundary conditions'), xlabel('x'), ylabel('y'), zlabel('z');
colorbar, colormap cool;
shading interp;
set(f,'EdgeColor',[0 0 0]);

微分方程式を離散化して解く

この例では、有限差分近似法を使用して、微分方程式を線形システムに離散化します。また、マルチグリッド前処理行列を使用して、反復ソルバーのパフォーマンスを向上させます。この例では、サポート関数 discretizePoissonEquation および multigridPreconditioner で離散化と前処理行列を使用します。その他の問題については、アプリケーションに適切な離散化と前処理行列を選択してください。

この例では、discretizePoissonEquation はポアソン方程式を 7 点ステンシル有限差分法で、粒度のレベルが異なる複数のグリッドに離散化します。この関数は、事前計算済みの三角分解、および粒度の低いレベルと高いレベルの間でマッピングする演算子を使用して、離散化のマルチグリッド構造を作成します。前処理行列はこのマルチグリッド情報を使用して、解を効率的に近似します。

この前処理行列の特別な点として、平滑化を適用し、一連の近似で発生するエラーを最小限に抑えていることがあげられます。平滑化ステップの数を定義します。ステップ数を多くすると、近似の精度は高くなりますが、計算量が多くなります。その後、微分方程式を離散化して、前処理行列を設定します。

numberOfSmootherSteps = 1;
[A,b,multigridData] = discretizePoissonEquation(numPoints,multigridLevels,numberOfSmootherSteps,u);

Level 0: The problem is of dimension 32768 with 223232 nonzeros.
Level 1: The problem is of dimension 4096 with 27136 nonzeros.
Level 2: The problem is of dimension 512 with 3200 nonzeros.

preconditioner = setupPreconditioner(multigridData);

前処理付き共役勾配法を使用して線形システムを解きます。

tol = 1e-12;
maxit = numel(A);
pcg(A,b,tol,maxit,preconditioner);

pcg converged at iteration 45 to a solution with relative residual 5.4e-13.

分散配列によるスケールアップ

メモリなどの計算リソースがさらに多く必要となる場合、分散配列を使用して、コードを変更せずにスケールアップできます。分散配列は、データを複数のワーカー間で分散します。また、マシンのクラスターの計算パフォーマンスとメモリを利用できます。

並列ワーカーのプールを起動します。既定で、parpool は既定のクラスターを使用します。既定のクラスタープロファイルを確認するには、MATLAB の [ホーム] タブの [環境] 領域で、[並列]、[既定のクラスターの選択] を選択します。

parpool;

Starting parallel pool (parpool) using the 'MyCluster' profile ...
Connected to the parallel pool (number of workers: 12).

関数 distributed を使用して、温度の変数 u をクラスターのワーカーのメモリ内に分散します。

distU = distributed(u);

コードは以前と同じものを使用できます。入力が分散配列の場合は離散化と前処理行関数により分散配列が作成されるため、変更は必要ありません。MATLAB 関数の多くは分散配列に対応しているため、メモリ内配列の場合と同様に作業ができます。

discretizePoissonEquation は分散データが含まれる構造体を返すことに注意してください。構造体内の分散データを分散して使用するには、spmd ブロック内で構造体を作成しなければなりません。また、これを使用する関数はすべて spmd ブロック内で呼び出さなければなりません。

spmd
    [A,b,multigridData] = discretizePoissonEquation(numPoints,multigridLevels,numberOfSmootherSteps,distU);
    preconditioner = setupPreconditioner(multigridData);
end

Analyzing and transferring files to the workers ...done.
Lab  1: 
  Level 0: The problem is of dimension 32768 with 223232 nonzeros.
  Level 1: The problem is of dimension 4096 with 27136 nonzeros.
  Level 2: The problem is of dimension 512 with 3200 nonzeros.

spmd ブロック内で pcg を使用し、線形システムを分散して解きます。

spmd
    x = pcg(A,b,tol,maxit,preconditioner);
end

Lab  1: 
  pcg converged at iteration 45 to a solution with relative residual 5.4e-13.

結果のプロット

ソルバーによる解は、メモリに収まるベクトルです。gather を使用して、ワーカーからクライアントにデータを送信します。データの形状を 3 次元配列に戻し、最終的な解を生成するために次元を並べ替えます。u の内側の部分をこの解に設定します。外側の部分 (境界) には境界条件の値が既に含まれています。

x3D = reshape(gather(x),numPoints.X,numPoints.Y,numPoints.Z);
u(2:end-1,2:end-1,2:end-1) = permute(x3D, [2, 1, 3]);

関数 slice を使用して解を可視化します。スライスを追加して、立方体内の温度をプロットします。解を調べるには、Rotate ツールを使用するか、スライスの位置を変化させます。

xSlices = [.5,1];
ySlices = [.5,1];
zSlices = [0,.5];
f = slice(X,Y,Z,u,xSlices,ySlices,zSlices,'nearest');
title('Heat distribution'), xlabel('x'), ylabel('y'), zlabel('z');
colorbar, colormap cool;
shading interp;

この例では numPoints のさまざまな値を試して、複数の離散化レベルをテストすることができます。大きな値を使用すると、分解能は高くなりますが、必要となるメモリも多くなります。また、multigridLevels の値が大きいと、前処理行列のメモリ効率が高くなります。ただし、multigridLevels の値が大きくなると、粒度が低くなってすべてのレベルで精度が下がるため、前処理行列の精度が低くなります。その結果、同じ水準の精度を実現するにはソルバーで多くの反復が必要となる場合があります。

前処理行列の定義

前処理付き共役勾配法で使用するマルチグリッド前処理行列を定義します。このような前処理行列では、粒度が異なる複数の離散化グリッドを使用して、線形方程式系の解を効率的に近似します。この例の前処理法は [2] に基づいており、以下の主な段階に従います。

ガウス - ザイデル近似法を使用して事前に平滑化する。
粗いレベルで残差解を計算する。
粗いレベルで再帰的に前処理を実行するか、最も粗いレベルの場合は直接解く。
粗いグリッドの解を使用して解を更新する。
ガウス - ザイデル近似法を使用して事後に平滑化する。

function x = multigridPreconditioner(mgData,r,level)

if(level < mgData(level).MaxLevel)
    x = zeros(size(r),'like',r);
    
    % Presmooth using Gauss-Seidel
    for i=1:mgData(level).NumberOfSmootherSteps
        x = mgData(level).Matrices.LD \ (-mgData(level).Matrices.U*x + r);
        x = mgData(level).Matrices.DU \ (-mgData(level).Matrices.L*x + r);
    end
    
    % Compute residual on a coarser level
    Axf = mgData(level).Matrices.A*x;
    rc = r(mgData(level).Fine2CoarseOperator)- Axf(mgData(level).Fine2CoarseOperator);
    
    % Recursive call until coarsest level is reached
    xc = multigridPreconditioner(mgData, rc, level+1);
    
    % Update solution with prolonged coarse grid solution
    x(mgData(level).Fine2CoarseOperator) = x(mgData(level).Fine2CoarseOperator)+xc;
    
    % Postsmooth using Gauss-Seidel
    for i = 1:mgData(level).NumberOfSmootherSteps
        x = mgData(level).Matrices.LD \ (-mgData(level).Matrices.U*x + r);
        x = mgData(level).Matrices.DU \ (-mgData(level).Matrices.L*x + r);
    end
else
    % Obtain exact solution on the coarsest level
    x = mgData(level).Matrices.A \ r;
end

end

マルチグリッドデータを取り前処理行列を入力データに適用する関数ハンドルを返す関数を作成します。この例では、この関数ハンドルは pcg への前処理行列入力です。spmd ブロック内に無名関数を定義することはできないため、この関数は必ず作成しなければなりません。

function preconditioner = setupPreconditioner(multigridData)

if ~isempty(multigridData)
    preconditioner = @(x,varargin) multigridPreconditioner(multigridData,x,1);
else
    preconditioner = [];
end

end

参考文献

[1] Dongarra, J., M. A. Heroux, and P. Luszczek. HPCG Benchmark: A New Metric for Ranking High Performance Computing Systems. Knoxville, TN: University of Tennessee, 2015.

[2] Elman, H. C., D. J. Silvester, and A. J. Wathen. Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics. Oxford, UK: Oxford University Press, 2005, Section 2.5.

参考

pcg | distributed | spmd | gather