読者です 読者をやめる 読者になる 読者になる

lochtext

http://us.battle.net/sc2/en/profile/2312833/1/lochtext/

移流方程式

久々に数値計算でもするかーと思い仮想環境を入れたら思ったよりうまく行ってうれしかった。
VertualBox+Lubuntu、エディタはvimを使ってみた。

何するかーってなって、今回は
Fortran 勉強会 wiki - itani-1 次元線形移流方程式
の先頭のやつをC++で書きなおした。べつにCでも良かったと思う。


#include <iostream>
#include <math.h>
#include <fstream>
#include <string>
using namespace std;

int main(int argc,char *argv[]){

 const int jmax = 100;       //  0 < x < 10   で計算 
 const int nmax = 500;       //  0 < t < 10 s で計算
 class p{
   public:
   double x;
   double u[nmax+1];
 };

 p point[jmax+1];
 
  const double pi = 3.14;
  const double c  = 1.0;
  const double dx = 0.1;
  const double dt = 0.02;
  double r  = c * (dt/dx);

//--- write parameters -----
  std::cout <<"dt="<<dt<<"dx="<<dx<<"r="<<r<<std::endl;

//--- initialize x -----
  point[0].x = 0.0;


//--- set coordinate -----
  for(int j = 1; j <= jmax; j++){
    point[j].x = point[j-1].x + dx;
    //cout << j << "  " << point[j].x << endl;
  }

//!--- set initial (t=1) condition -----
//! sin 波と矩形波の二種類用意した.
//! 与えた値によってスイッチする.

  cout <<
  "*****************************************"<<endl<<
  "*       choose initial condition.       *"<<endl<<
  "*  1 (= sin wave) or 2 (= square wave)  *"<<endl<<
  "*****************************************"<<endl;

  int iniwave;
  cin >> iniwave;

  ofstream fout("dat/000.dat");
  if(iniwave==1){

    //! sin 波 : 0 < x < pi で sin(x)
    for(int j = 0; j <= jmax; j++){
      if(point[j].x < pi){
        point[j].u[0] = sin(point[j].x);
      }
      else{
        point[j].u[0] = 0.0;
      }
      fout << point[j].x << "  " << point[j].u[0] << endl; 
    }
  }
  else{
    //! 矩形波 : u = 0 ( x < 1, 2 < x ), u = 1.0 ( 1 <= x <= 2  )
    for(int j = 0; j <= jmax; j++){
      if(point[j].x>=1.0 && point[j].x<=2.0){
        point[j].u[0] = 1.0;
      }
      else{ 
        point[j].u[0] = 0.0;
      }
      fout << point[j].x << "  " << point[j].u << endl; 
    }
    fout.close();
  }


//--- calculate each value -----


  for (int n = 1; n <= nmax ; n++){
    char name[30];
    sprintf(name,"dat/%03d.dat",n);
    ofstream fout(name);
        
    point[0].u[n]= point[jmax].u[n-1];
    fout << point[0].x << "  " << point[0].u[n] << endl;

    for (int j = 1 ; j <= jmax ; j++){    
      point[j].u[n] = (1.0-r) * point[j].u[n-1] + r * point[j-1].u[n-1];
      fout << point[j].x << "  " << point[j].u[n] << endl;
    }
    fout.close();
  }

return 0;

}


載ってたシェルスクリプトをちょこっと変えて、Gnuplot->PNG->ImageMagickでconvert、GIFアニメ作成までできた。シェルスクリプトは割愛。

シェルスクリプトGnuplotC++と、やらないとどんどん忘れていくのでたまにやるようにしたい。あと、流体解析の基本的なところは卒業までに身につけたい。FortranC++で書き換えるだけで勉強した感じになれてよい。

次回はPythonでやるか、Eigenを使ってCrank-Nicolson法でも作ろうかな。