Eigenを用いた2次元移流拡散方程式の数値計算

二次元空間での移流拡散方程式の数値計算です。拡散する物理量は濃度とか温度だと思ってください。

離散化手法として、時間離散化にクランク＝ニコルソン法、空間離散化に１次精度の風上差分法を用いています。

二次元での陰解法なので、TDMAで線順法（c.f.パタンカー）でもいいのですが、Eigenの疎行列ライブラリが使いたかったので一次方程式にして解きました。あんまりオブジェクト指向してませんがそこは適当です。

連立方程式のソルバーは適当に選んでます。このへんの知識をもう少し入れたいですね。

diffusion.cpp

#include <iostream>
#include <math.h>
#include <fstream>
#include <string>
#include <Eigen/Sparse>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

typedef SparseMatrix<double> SpMat;

int main(){

  const int imax = 100;//x方向の分割数
  const int jmax = 100;//y方向の分割数
  const int nmax = 100;//時間分割数
  
  const double xmax = 10.0;//x方向の距離
  const double ymax = 10.0;//y方向の距離
  const double Time = 10.0;//z方向の距離

  class position{//2次元空間の升目ごとにインスタンスを作成
    public:
    int i;//number(x)
    int j;//number(y)
    double x;//position
    double y;//position
    double *u;//拡散する物理量へのポインタ。
    double *uold;
    position *E, *W, *N, *S;//二次元空間で上下左右のマスへのポインタ
  };

  position p[imax+2][jmax+2];//升目の数だけインスタンスを作成。境界条件はi=0 or i=imax+1でu=0。
  double u[imax+2][jmax+2];//拡散する物理量。
  double uold[imax+2][jmax+2];//時間発展用。

  const double pi = 3.14;//円周率
  const double C  = 0.5;//移流係数
  const double D  = 0.1;//拡散係数
  double dx = xmax/double(imax);
  double dy = ymax/double(jmax);
  double dt = Time/double(nmax);
  double c = dt * C /dx;
  double d = dt * D /dx /dx;

//set coordinate
  for(int i = 0; i <= imax+1; i++){
    for(int j = 0; j <= jmax+1; j++){
      if(i==0){
        p[i][j].x = 0.0;
      }else{
        p[i][j].x = p[i-1][j].x+dx;
      }
      
      if(j==0){
        p[i][j].y = 0.0;
      }else{
        p[i][j].y = p[i][j-1].y+dy;
      }
      p[i][j].i=i;
      p[i][j].j=j;
      p[i][j].E = (i!=imax+1)? &p[i+1][j]: &p[i][j];
      p[i][j].W = (i!=0)?      &p[i-1][j]: &p[i][j];;
      p[i][j].N = (j!=jmax+1)? &p[i][j+1]: &p[i][j];;
      p[i][j].S = (j!=0)?      &p[i][j-1]: &p[i][j];;
      p[i][j].u = &u[i][j];
      p[i][j].uold = &uold[i][j];
    }
  }

  ofstream fout("dat/000.dat");


//初期条件として2次元のsin波をセット
  for(int i = 0; i <= imax+1; i++){
  for(int j = 0; j <= jmax+1; j++){
    double center = 1.0;
    double d = sqrt((center-p[i][j].x)*(center-p[i][j].x)+(center-p[i][j].y)*(center-p[i][j].y));
    if(i==0||i==imax||j==0||j==jmax){
      uold[i][j] = 0.0;
    }else if(d < pi/2 ){
      uold[i][j] =sin(pi/2-d);
    }else{
      uold[i][j] = 0.0;
    }
    if(i%2==0&&j%2==0){fout << p[i][j].x << "  "<< p[i][j].y << "  "  << uold[i][j] << endl;}
  }}

  fout.close();
  
  cout << "000.end" <<endl;


//implicit method

  int mmax = (imax+2)*(jmax+2);//行列の大きさ
  SparseMatrix<double> A(mmax,mmax),B(mmax,mmax);
  MatrixXd b(mmax,1),xb(mmax,1),x(mmax,1);
  BiCGSTAB<SpMat> solver;

//二次元の移流拡散方程式を陰開放で離散化すると、
//最終的に A*u(n+1)=B*u(n)という
//行列を用いた形にできます。以下ではAおよびBを作成しています。
//問題によっては（たとえば対流項が非線形なNS方程式など）
//時間発展ごとにA,Bを再作成する
//必要がありますが、今回は一定でOKです。

  //insert to A
  vector<Triplet<double> > t1;
  for(int i=0;i<=imax+1;i++){
  for(int j=0;j<=jmax+1;j++){
    t1.push_back(Triplet<double>((jmax+2)*i+j,(jmax+2)*i+j,1-c+2*d));
    t1.push_back(Triplet<double>((jmax+2)*i+j,(jmax+2)*p[i][j].E->i+p[i][j].E->j,c/2-d/2));
    t1.push_back(Triplet<double>((jmax+2)*i+j,(jmax+2)*p[i][j].W->i+p[i][j].W->j,-d/2));
    t1.push_back(Triplet<double>((jmax+2)*i+j,(jmax+2)*p[i][j].N->i+p[i][j].N->j,c/2-d/2));
    t1.push_back(Triplet<double>((jmax+2)*i+j,(jmax+2)*p[i][j].S->i+p[i][j].S->j,-d/2));
  }}
  A.setFromTriplets(t1.begin(),t1.end());
  solver.compute(A);

  //insert to B
  vector<Triplet<double> > t2;
  for(int i=0;i<=imax+1;i++){
  for(int j=0;j<=jmax+1;j++){
    t2.push_back(Triplet<double>((jmax+2)*i+j,(jmax+2)*i+j,1+c-2*d));
    t2.push_back(Triplet<double>((jmax+2)*i+j,(jmax+2)*p[i][j].E->i+p[i][j].E->j,-c/2+d/2));
    t2.push_back(Triplet<double>((jmax+2)*i+j,(jmax+2)*p[i][j].W->i+p[i][j].W->j,d/2));
    t2.push_back(Triplet<double>((jmax+2)*i+j,(jmax+2)*p[i][j].N->i+p[i][j].N->j,-c/2+d/2));
    t2.push_back(Triplet<double>((jmax+2)*i+j,(jmax+2)*p[i][j].S->i+p[i][j].S->j,d/2));
  }}
  B.setFromTriplets(t2.begin(),t2.end());

  cout << "insert to A&B.end" <<endl;

//main loop

  for (int n = 1; n <= nmax ; n++){
     
    char name[30];
    sprintf(name,"dat/%03d.dat",n);
    ofstream fout(name);
    
    if(n%10 == 0){
      cout << "this is step" << n << endl;
    }
    
    b = Map<MatrixXd> (&uold[0][0],mmax,1);
    xb = B*b;
    x = solver.solveWithGuess(xb,b);
    
    Map<MatrixXd>(&u[0][0],mmax,1) = x;

    for (int i = 1;i<=imax;i++){
    for (int j = 1;j<=jmax;j++){
      uold[i][j] = u[i][j];
      
      if(i%2==0&&j%2==0){
        fout << p[i][j].x << "  " << p[i][j].y << "  " << u[i][j] << endl;
      }
    }}
    fout.close();

  }

}

C++ではなくPythonですが、数値流体力学の基礎を学べるコンテンツがあったので紹介しておきます。陽解法でキャビティ流れを解くのが最終目標のようです。