Fortranで数値計算の結果をデータベースで管理する

概要

以前に『Fortranでもデータベースしたい！』という題で記事を投稿しました。本稿ではそれに続いて、数値計算プログラムのデータ管理をLibpq-FortranとPostgreSQLで行う方法について例示したいと思います。

Libpq-Fortranとは

データベースサーバーPostgreSQLへのFortranインターフェースです。

Fortran2003以降のC相互運用機能を利用してPostgreSQLのCライブラリlibpqの関数をFortranから呼び出せるようにしました。libpqのFortranラッパーとも言えるでしょう。

これを利用するとFortranプログラムからPostgreSQLのサーバーに接続してクエリを実行し、データを読み書きすることができます。

GitHubで公開中です

数値計算のデータ管理

Fortranプログラムを扱う人のその目的は、多くの場合、数値計算プログラムを書いて結果を評価することでしょう。数値シミュレーションを例に、このプロセスをデータに注目して分けると以下のようになります（各サブプロセスで扱うデータファイルを括弧内に書いています）。

1. ソースコードを書く（ソースファイル）
2. シミュレーションを実行する（入力ファイルと出力ファイル）
3. 結果の図をプロット（描画スクリプトと画像ファイル）
4. 1に戻る

さて、これらのファイルはどのように取り扱えば良いでしょうか？

1番目のソースコードや3番目のスクリプトの管理には、GitやSVNなどのバージョン管理システムを使用することができます。

しかし一方で、2番のシミュレーションの初期値（入力ファイル）や計算結果（出力ファイル）を管理するためのベストプラクティスは自明ではありません。単純にファイルシステムで管理すると仮定すれば、以下のような方式が考えられるでしょう。

• 識別できるファイル名をつける
  • 例：hogefuga_202312022400_A.dat
    • 実行した日時などを入れる。
    • できるだけキーワードを入れて長い名前にする。
• ディレクトリによるカテゴリ分けをする。
  • 例：result/hoge-scheme/fugapiyo_20231202.dat

しかしこれらの方法も、並列処理をすると衝突してしまったり、ディレクトリパスが長くなってしまったりして、データを取り出すのに手間取ることになります。また、そもそも長い名前を考えるのは多少の面倒でもあります。

このような手間は、データ数が小さいうちは手動で間に合いますが、数が膨大になると管理しきれない場面が出てくることでしょう。

入出力データのデーターベース管理

前のセクションで、数値計算の結果を管理するのは悩みの種になることについて述べました。このセクションではデータベース管理システムの一つであるPostgreSQLを利用して、それらのデータを管理する方法について詳しく見ていきます。

Libpq-Fortranを使用すると、Fortranプログラムからデータベースサーバーにアクセスすることが可能になります。以下では実際にPostgreSQLデータベースにアクセスして、どのように数値計算の結果を登録するかについて例示します。

例：一次元移流方程式の数値シミュレーション

一次元移流方程式の数値シミュレーションを例にとり、実際にパラメーターを色々変えて実行し、結果をデータベースに登録してみましょう。

一次元移流方程式

$$\frac{\partial u}{\partial t} + \frac{\partial f}{\partial x} = 0, f= cu,$$

を陽的オイラー法と、Lax-FriedrichスキームおよびLax-Wendroffスキームでそれぞれ離散化すると、次のような差分方程式が得られます。

$$
u^{n+1}_{j}= u^{n}_{j} - \frac{\Delta t}{\Delta x}\left(\tilde f^{n+\frac{1}{2}}_{j+\frac{1}{2}}-\tilde f^{n-\frac{1}{2}}_{j-\frac{1}{2}}\right)
$$

ここで$\tilde f$は数値流速であり、各スキームで以下のように取ります。

Lax-Friedrichスキーム：

$$
\tilde f^{n+\frac{1}{2}}_{j+\frac{1}{2}}=\frac{1}{2}\left[ \left(f^n_{j+1}+f^n_j\right)-\frac{\Delta x}{\Delta t}\left(u^n_{j+1}-u^n_j\right)\right],
$$

Lax-Wendroffスキーム：

$$
\tilde f^{n+\frac{1}{2}}_{j+\frac{1}{2}} = \frac{1}{2}\left[\left(1 - c\frac{\Delta t}{\Delta x}\right)f^n_{j+1}+ \left(1+ c\frac{\Delta t}{\Delta x}\right)f^n_j\right],
$$

なお、ここで$$f^n_{j} = f(u^n_j) = cu^n_j$$です。

また、$c$は特性速度、$\Delta t$は時間ステップ幅、$\Delta x$は空間刻み幅と呼ばれる量になり、この方程式のシミュレーション1回を特徴づけるパラメーターになります。

ユーザーとデータベースを作る

以下のFortranプログラムでは、PostgreSQLサーバーに接続して、2つのクエリを実行しています。

program create_database
   use, intrinsic :: iso_c_binding
   use :: libpq
   implicit none
   
   type(c_ptr) :: conn, res
   character(:), allocatable :: conninfo, query
   
   conninfo = "host=localhost user=postgres"  ! 接続情報
   !（必要なら管理者ユーザーpostgresのパスワードをpassword=以下に記述する）
   
   conn = PQconnectdb(conninfo)               ! 接続の実行

   ! エラーハンドリング
   if (PQstatus(conn) /= CONNECTION_OK) then
      print *, PQerrorMessage(conn)
      error stop
   end if

   ! ユーザー作成のクエリ
   query = "create role shinobu with createdb login password 'hogehoge';"
   res = PQexec(conn, query)    ! クエリ実行

   ! エラーハンドリング
   if (PQresultStatus(res) /= PGRES_COMMAND_OK ) then
      print *, PQresultErrorMessage(res)
      call PQfinish(conn)
      error stop
   end if
   
   ! 結果オブジェクトの解放
   call PQclear(res)

   ! データベース作成のクエリ
   query = "create database simulation with owner shinobu;"
   res = PQexec(conn, query)    ! クエリの実行

   ! エラーハンドリング
   if (PQresultStatus(res) /= PGRES_COMMAND_OK ) then
      print *, PQresultErrorMessage(res)
   end if

   call PQclear(res)    ! 結果オブジェクトの解放
   call PQfinish(conn)  ! 接続オブジェクトの解放

end program create_database

データベースへ接続

最初にconninfo変数に接続情報（ホスト名localhost、ユーザー名postgres）を記述し、関数PQconnectdbを呼び出して、戻り値をtype(c_ptr)の変数connに代入しています。

その後のif文はPQstatus(conn) /= CONNECTION_OKという条件文を評価して、エラーハンドリングを行っています。

ユーザー作成のクエリ

次にユーザー作成のクエリの部分では、create role shinobu with createdb login password 'hogehoge';という文字列をquery変数に格納した後、PQexec関数を呼び出して、その結果をtype(c_ptr)型変数resに代入しています。

ここでquery変数の中身にはデータベースサーバーに対するコマンドを含んでいて、以下のような意味を持ちます。

• create role shinobu：shinobuという名前でロール（ユーザー）を作成するコマンド
• with createdb：当該ロールにデータベースを作成する権限を付与する。
• login password 'hogehoge'：当該ロールのログインパスワードをhogehogeに設定する。

ここで、ユーザー名とパスワードは例示したものではなく各自で置き換えてください。

PQexec関数は、第一引数にPQconnectdbで返された接続オブジェクトへのポインタを、第二引数に上記のクエリを取り、結果オブジェクトへのポインタtype(c_ptr)型の値を返します。ここではtype(c_ptr)型の変数resに代入しています。

コマンドが成功した場合はPQresultStatus(res)がPGRES_COMMAND_OKの値を返します。失敗した場合にはこれ以外の値を返し関数PQresultErrorMessage(res)がエラーメッセージを返します。

成功した場合は、もうresオブジェクトは不要なのでPQclear(res)関数を呼び出して、このオブジェクトを解放します。

データベース作成のクエリ

次にデータベース作成のクエリの文字列create database simulation with owner shinobu;をquery変数に格納しています。

ここでクエリに含まれるコマンドの意味は次のようなものです：

• create database simulation：simulationという名のデータベース（テーブルの集合の名前）を作成する。
• with owner shinobu：当該データベース（simulation）の所有者をshinobuにする。

その後、この文字列をPQexec関数に渡して呼び出しています。エラーハンドリングについては上と同様です。

成功した場合、結果のオブジェクトresと接続オブジェクトconnを解放して、このプログラムは終了です。

この時点で、ターミナルからpsql -h localhost -U shinobu -d simulaitonを実行して、ここで作成したデータベースsimulationに接続することができるはずです。

テーブルを作る

次に、実際にシミュレーションの情報を書き込むためのテーブルを用意します。

conninfoに記述する接続情報には、前のセクションで作成したロール・データベース・パスワードを使用します。

program create_table
   use :: libpq
   use, intrinsic :: iso_c_binding
   implicit none
   
   type(c_ptr) :: conn, res
   character(:), allocatable :: conninfo, query
   integer :: i

   conninfo = "host=localhost user=shinobu dbname=simulation password=hogehoge"
   conn = PQconnectdb(conninfo)
   if (PQstatus(conn) /= 0) then
      print *, PQerrorMessage(conn)
      error stop
   end if

   query = "create table advection &
   & (sim_id      decimal(4,0)      not null, &
   &  start_date  date              ,         &
   &  sim_name    varchar(100)      not null, &
   &  cs          double precision  not null, &
   &  dt          double precision  not null, &
   &  dx          double precision  not null, & 
   &  courant     double precision  ,         &
   &  length      double precision  ,         &
   &  nx          integer           not null, &
   &  nstep       integer           not null, &   
   &  dat_file    varchar(256)      not null, &
   &  primary key(sim_id));"
   
   res = PQexec(conn, query)
   if (PQresultStatus(res) /= PGRES_COMMAND_OK) then
      print *, PQresultErrorMessage(res)
   end if 

   call PQclear(res)
   call PQfinish(conn)

end program create_table

変数queryには最初にcreate table advectionとあり、ここでadvectionという名前のテーブルを定義することを宣言しています。その後、括弧の中に属性名・型・オプションを1つの属性（テーブルの列）としてカンマ区切りで定義を並べています。各属性の記述の意味は次の通りです。

• sim_id: 実験ID、整数型（4桁を小数点以下0桁精度）、非NULL制約¹あり
• start_date：実行した日付、日付型
• sim_name：シミュレーションの名前、最大100桁の可変長文字列型、非NULL制約あり
• cs：特性速度の数値、倍精度実数型、非NULL制約あり
• dt：時間ステップ幅、倍精度実数型、非NULL制約あり
• dx：空間刻み幅、倍精度実数型、非NULL制約あり
• courant：クーラン数=$c\Delta t/\Delta x$、倍精度実数型
• length：計算領域長さ、倍精度実数型
• nx：計算領域の分割数、整数型、非NULL制約あり
• nstep：計算した時間ステップの総数、非NULL制約あり
• dat_file：データファイルのパス名、最大256桁の可変長文字列型、非NULL制約あり

最後にこのテーブルでは一意にレコードを特定できる属性として属性sim_idを選ぶことにするので、一意性制約primary key(sim_id)を記述しています。

これらを一つの変数queryに含ませた状態で、PQexecを実行します。

このプログラムを実行した後に、ターミナルからpsql -h localhost -U shinobu -d simulationを実行した上で、psqlのプロンプトからselect * from advection;を問い合わせて、以下のような属性名の一覧を得られれば成功です。

$ psql -h localhost -U shinobu -d simulation
psql (15.4)
Type "help" for help.

simulation=> select * from advection;
 sim_id | start_date |   sim_name    | cs |   dt   |  dx   | courant | l | nx  | nstep |      dat_file
--------+------------+---------------+----+--------+-------+---------+---+-----+-------+--------------------

以下のセクションではこのテーブルに対してデータの登録や取り出しを試していきます。

シミュレーションを実行する。

データを保持しておくデータベースとテーブルが用意できたので、シミュレーションを実行していきます。

ここで実行するプログラムは以下の要素から構成されています。

• 主プログラム：1個の実験をパラメータを変更して全部で32個実行するループを持つ。
• 内部サブルーチン
  • main_loopサブルーチン：引数にシミュレーションを特徴づけるパラメータを受け取り、1個のシミュレーションを実行する。
  • data_outputサブルーチン：出力ファイルの第1列にx、第2列にuの値を書き出し、区切りを2つの空行とするデータセットを書き込む。
  • data_registrationサブルーチン：main_loopの後に、実験のパラメータと出力ファイルへのパスをデータベースに登録する

時間積分するメインのループ

数値計算のループはmain_loopというサブルーチンに含まれているものとします。このサブルーチンの引数に、1回の実験を特徴づけるパラメータを渡して、変数filenameにはファイル名を作成して格納し、その名前の.datファイルに結果を書き込みます。

subroutine main_loop(cs, dt, dx, nstep, filename, offset_k, k, prefix)
   implicit none
   real(real64), intent(in) :: cs, dt, dx ! cs is short for characteristic speed
   integer(int32), intent(in) :: k
   character(:), allocatable, intent(inout) :: filename
   integer :: offset_k
   character(*), intent(in) :: prefix
   integer, intent(out) :: nstep

   character(4) :: code
   integer :: nth_out
   real(real64) :: tout, time
   real(real64) :: Courant

   integer :: i
   integer :: uni

   real(real64), allocatable :: x(:)
   real(real64), allocatable :: u(:), f(:), u_new(:), f0(:)

   write(code, '(1i4.4)') k + (offset_k - 1)

   Courant = cs*dt/dx

   write(error_unit, '(a, i0, a)') "### START ", k, "-th Simulation ###"
   write(error_unit, '(a, f10.3)') " Courant number: ", Courant

   ! Initialize
   uni = 100
   nth_out = 0
   nstep = 1
   time = 0d0
   tout = 0
   allocate(x(0:nx), source=0d0)
   allocate(u(0:nx), source=0d0)
   allocate(u_new(0:nx), source=0d0)
   allocate(f0(0:nx), source=0d0)
   allocate(f(0:nx), source=0d0)
   
   ! 空間のグリッドをセットする
   x(0) = 0d0
   do i = 1, nx
      x(i) = x(0) + dx*dble(i)
   end do

   ! 初期条件をセットする
   do i = 0, nx
      if (x(i) < 0.1d0) then
         u(i) = 1d0
      else
         u(i) = 0d0
      end if
   end do

   filename = 'result/'//trim(adjustl(prefix))//code//'.dat'
   
   open(uni, file=filename, form='formatted', status='replace')
   call data_output(uni, x, u, tout, time, 0, nth_out)

   ! 時間積分
   do
      nstep = nstep + 1
      time = time + dt

      f0(:) = cs*u(:) ! 数値流束

      ! Lax-Friedrich スキーム
      if (trim(prefix) == 'LF_') then
         do i = 1, nx-1
            f(i) = 0.5d0*(f0(i+1)+f0(i)) - 0.5d0*dx/dt *(u(i+1) - u(i))
         end do          
      end if

      ! Lax-Wendroff スキーム
      if (trim(prefix) == 'LW_') then
         do i = 1, nx-1
            f(i) = 0.5d0*( (1d0-Courant)*f0(i+1) + (1d0+Courant)*f0(i) )
         end do
      end if

      ! 陽的オイラー法
      do i = 1, nx-1 
         u_new(i) = u(i) - Courant*(f(i) - f(i-1))
      end do

      ! uのアップデート
      u(2:nx-1) = u_new(2:nx-1)

      ! 境界条件
      u(1) = u(2)
      u(nx) = u(nx-1)

      if (time >= tout) call data_output(uni, x, u, tout, time, nstep, nth_out)
      if (time >  tend) exit

   end do

   close(uni)

   write(error_unit, fmtStop) nstep, time
   write(error_unit, '(a)') "### Normal STOP ###"

end subroutine main_loop

main_loopサブルーチンに含まれるdoループ内部に、セクションの冒頭で式で示したLax-FriedrichスキームとLax-Wendroffスキームの数値流束の計算を引数prefixの値に応じて切り替えて実行するようにしています。

ファイル出力

データ書き込み用のサブルーチンは以下の通りです。

subroutine data_output (uni, x, u, tout, time, nstep, nth_out)
   implicit none
   integer, intent(in) :: uni
   real(real64), intent(in) :: x(:), u(:)
   real(real64), intent(inout) :: tout
   real(real64), intent(in) :: time
   integer, intent(in) :: nstep
   integer, intent(inout) :: nth_out 
   integer :: i

   do i = 1, nx
      write(uni, fmtOut) x(i), u(i)
   end do
   write(uni, '(/)')  !空行を2つ書き込む(gnuplot用のデータセットの区切り) 

   tout = tout + tout_interval
   
   write(error_unit, fmtErrOut) nstep, time, nth_out
   nth_out = nth_out + 1
end subroutine data_output

データを登録する

データベースにデータを登録するためのサブルーチンは次の通りです。ここでのクエリはinsert文と呼ばれるもので、テーブルにレコードを挿入します。

subroutine data_registration(conn, simu_name, filename, offset_k, k, cs, dt, dx, nstep, date)
   use :: iso_c_binding
   use :: libpq
   implicit none
   type(c_ptr), intent(in) :: conn
   character(*), intent(in) :: simu_name
   character(*), intent(in) :: filename
   integer, intent(in) :: k, offset_k
   real(real64), intent(in) :: cs, dt, dx
   integer, intent(in) :: nstep
   character(*), intent(in) :: date

   character(:), allocatable :: query
   character(256) :: values(11)

   type(c_ptr) :: res
   real(real64) :: Courant
   character(4) :: code
   integer :: i

   Courant = cs*dt/dx

   ! 文字列配列に値を格納する
   write(code, '(1i4.4)') k + (offset_k-1)
   values(1) = code
   values(2) = trim(adjustl(date))
   values(3) = trim(simu_name)
   write(values(4),  '(f16.8)') cs
   write(values(5),  '(f16.8)') dt
   write(values(6),  '(f16.8)') dx
   write(values(7),  '(f16.8)') Courant
   write(values(8),  '(f16.8)') L 
   write(values(9),  '(i0)') nx
   write(values(10), '(i0)') nstep
   values(11) = trim(adjustl(filename)) 

   do i = 1, 9
      values(i) = trim(adjustl(values(i)))
   end do

   ! プレースホルダーには文字列配列valuesの要素の値が展開される
   query = 'insert into advection &
   &        (sim_id, start_date, sim_name, cs, dt, dx, courant, l, nx, nstep, dat_file) &
   & values ($1,     $2,         $3,       $4, $5, $6, $7,     $8, $9, $10,   $11)'

   res = PQexecParams(conn, query, 11, [(0, i=1,11)], values(:))

   if (PQresultStatus(res) /= PGRES_COMMAND_OK) then
      print *, PQresultErrorMessage(res)
   end if

   call PQclear(res)

end subroutine data_registration

データベースへの登録（データの挿入）を行うには、PQexecParams関数を使用してクエリを実行します。このPQexecParams関数は、クエリ文字列にはドル記号と数字からなるプレースホルダー（コード例では$1 から$11）を記述して、そこの値は文字列の配列として別に渡すことができます。PQexecParams関数に渡す引数の内容は次の通りです。

PQexecParams関数	引数の中身
第1引数	接続オブジェクトへのポインタ（conn）
第2引数	パラメータのプレースホルダーを含むクエリ文字列（query）
第3引数	第2引数に含まれるパラメータの個数（整数型）
第4引数	パラメータの個数だけの長さを持つゼロ値配列 （ここでのゼロはパラメータの型の解釈をサーバーに委ねるの意味）
第5引数	パラメータの個数だけの長さを持つ文字列型の配列

data_registrationサブルーチンの24行目から34行目の間に、第5引数にわたす予定の文字列の配列value(1:11)にデータの値を書き込んでいます。

メインプログラムを実行する

ここまでで必要な材料はすべて揃いました。メインプログラムを書いてシミュレーションを実行します。メインプログラムの全体はGitHubのこちらのリンクから入手することができます。

program main
   use, intrinsic :: iso_fortran_env
   use, intrinsic :: iso_c_binding
   use :: libpq

   implicit none

   real(real64), parameter :: tend = 0.6d0
   real(real64), parameter :: tout_interval = 0.1d0

   integer, parameter :: NX = 200     ! 分割格子の数
   real(real64), parameter :: L = 1d0 ! 計算領域長さ

   character(*), parameter  :: fmtOut   = '(e10.3, 1x, e10.3, 1x, i0)'
   character(*), parameter  :: fmtErrOut= "('  output  ', 'step=', i8, ' time=', e10.3, ' nth_out =', i3)"
   character(*), parameter  :: fmtStop  = "('  stop    ', 'step=', i8, ' time=', e10.3)"
   
   character(:), allocatable :: filename
   integer :: k

   real(real64) :: cs, dt, dx
   integer :: nstep, offset_k

   type(c_ptr) :: conn
   character(:), allocatable :: conninfo
   character(10) :: date = '2023-12-02'
   
   conninfo = 'host=localhost dbname=simulation user=shinobu'
   conn = PQconnectdb(conninfo)

   if (PQstatus(conn) /= 0) then
      print *, PQerrorMessage(conn)
      error stop 
   end if

   ! Lax-Friedrichスキームの実験を16件実行する（id = 1:16）。
   offset_k = 1
   do k = 1, 16
      cs = 1d0
      dt = 5d-4  + 3d-4*dble(k-1)
      dx = L/dble(nx)
      call main_loop(cs, dt, dx, nstep, filename, offset_k, k, 'LF_')
      call data_registration(conn, 'Lax-Friedrich', filename, offset_k, k, cs, dt, dx, nstep, date)
   end do
   
   ! Lax-Wendroffスキームの実験を16件実行する（id = 17:32）。
   offset_k = 17
   do k = 1, 16
   	cs = 1d0
   	dt = 5d-4  + 3d-4*dble(k-1)
   	dx = L/dble(nx)
   	call main_loop(cs, dt, dx, nstep, filename, offset_k, k, 'LW_')
      call data_registration(conn, 'Lax-Wendroff', filename, offset_k, k, cs, dt, dx, nstep, date)
   end do
   
contains
　...
end program main

データベースに登録されているのを確認する

ターミナルからpsql -h localhost -U shinobu -d simulationを実行してpsqlを起動し、select * from advection;を実行すると32行のレコードが挿入されていることが確認できます。

以上で、数値計算の結果のファイルをデータベースサーバーに登録することができました。次のセクションでは登録されたデータを取り出してFortranから図を描画する方法を述べたいと思います。

データベースからデータを取り出して図にプロットする

PostgreSQLにアクセスできるプログラミング言語ならFortranでなくとも、セクションタイトルのことをできますが、ここではLibpq-Fortranでデータをサーバーから読み込んでgnuplotで図を描画する方法について説明したいと思います²。

ここでの方針は、simulationデータベースのadvectionテーブルから一部のレコードを抽出して、パス名だけを出力し、文字列型変数にスペース区切りで連結してプロットスクリプトのコマンドライン引数として渡します。

以下がplot.f90の例です。

program plot
   use :: libpq
   use, intrinsic :: iso_c_binding
   implicit none
   
   character(:), allocatable :: conninfo, query, cmd
   type(c_ptr) :: conn, res
   integer :: i

   conninfo = "host=localhost dbname=simulation user=shinobu"
   conn = PQconnectdb(conninfo)
   if (PQstatus(conn) /= 0) then
      print *, PQerrorMessage(conn)
      error stop
   end if

   ! 実行するコマンドの先頭を記述する
	cmd = "gnuplot -c plot.gp "

   query = "select dat_file from advection where sim_name = 'Lax-Wendroff' and courant > 0.55;"
   res = PQexec(conn, query)
   
   if (PQresultStatus(res) /= PGRES_COMMAND_OK) then
      print *, PQresultErrorMessage(res)
      call PQfinish(conn)
      error stop
   end if 

   do i = 0, PQntuples(res) - 1
      ! クエリで取り出したデータをcmdに追記していく
      cmd = cmd // ' '// PQgetvalue(res, i, 0)
   end do

   ! コマンドを実行していく
	print *, cmd
   call system(cmd)

   call PQclear(res)
   call PQfinish(conn)
   
end program plot

ここで、"select dat_file from advection where sim_name = 'Lax-Wendroff' and courant > 0.55;"というクエリは、以下のような意味を持ちます。

queryの一部分	意味
from advection	advectionテーブルから、
where sim_name ='Lax-Wendroff' and courant > 0.55	属性sim_nameがLax-Wendroffであり、かつcourantが0.55より大きいレコードを抽出し、
select dat_file	dat_file属性の値（すなわちパス名）を出力する。

このクエリの結果はc_ptr型変数resの示すオブジェクトから、libpqの関数を使用して取り出します。PQexecの後のdoループでは、0からPQntuples(res) - 1までのループを指示しています。そしてそのループの中で、PQgetvalue(res, i, 0)で第i番目のレコードの第0列のデータを取り出し、cmd変数に連結しています。

その後、systemサブルーチンを呼び出し、cmdに含まれるコマンドを実行しています。

実際に実行されるコマンドは次のようになるでしょう。

gnuplot -c plot.gp  result/LW_0024.dat result/LW_0025.dat result/LW_0026.dat result/LW_0027.dat result/LW_0028.dat result/LW_0029.dat result/LW_0030.dat result/LW_0031.dat result/LW_0032.dat

ここで、plot.gpは次のようなGnuplotスクリプトになります。

set term pngcairo size 1600,1200
set termoption noenhanced
set output "plot9.png"

set xlabel "x"
set ylabel "u"
set xrange [0:1]
set yrange [-0.5:2.0]

set multiplot layout 3,3

do for [i=1:9] {
   path = ARGV[i]
   set title path

   unset xlabel
   unset ylabel
   if ( i % 3 == 1) {
      set ylabel "u"
   }
   if ( i > 6 ) {
      set xlabel "x"
   }

   plot  path index 0 using 1:2 with lines title 't = 0.0', \
      path index 1 using 1:2 with lines title 't = 0.1', \
      path index 2 using 1:2 with lines title 't = 0.2', \
      path index 3 using 1:2 with lines title 't = 0.3', \
      path index 4 using 1:2 with lines title 't = 0.4', \
      path index 5 using 1:2 with lines title 't = 0.5', \
      path index 6 using 1:2 with lines title 't = 0.6
}

unset multiplot

このプログラム（plot.f90）を実行すると次のような画像を出力することができます。

以上で、Libpq-FortranでデータをPostgreSQLサーバーから読み込んで、gnuplotスクリプトで図を描画することができました。

おわりに

ライブラリLibpq-Fortranを使用して、PostgreSQLで数値シミュレーションのデータを管理し、そのデータをサーバーから取り出して使用する方法について述べました。

この方法は同じプログラムをパラメーターを変えて大量に実行するような場面で活躍すると思われます。筆者はある差分法スキームの安定性を評価するためにパラメーターの全探索をしたときに（手で管理するのを諦めて）データベースを活用する方法を選択したことがありました。そのときはLibpq-Fortranがなかったので別の方法（Pythonだったかな？）で行っていましたが、今では本稿で説明した通りFortranプログラムから直接登録することができるようになりました。

今回扱った例では、1個のシミュレーションが1つのレコードに対応していたので比較的簡単に管理することができました。一方、1つのシミュレーションに対して複数のファイルが出力されるような場合を想定すると、データ管理の難易度はやや上がります。これについてはまたいずれ書きたいと思います。

コードのテスト環境

• Windows 10, GFortran MinGw-W64 13.1.0, Gnuplot 6.0-rc1
• Gentoo Linux, GFortran 13.2.1, Gnuplot 5.4.6

1. この属性を空（NULL）にしてレコードを登録することができないという意味。 ↩

2. 描画処理もFortranで行いたい場合には、PLPlotなどのライブラリを使う必要があるでしょう。 ↩