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

## 概要

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

## Libpq-Fortranとは

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

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

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

[GitHubで公開中です](https://github.com/ShinobuAmasaki/libpq-fortran)

## 数値計算のデータ管理

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

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

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

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

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

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

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

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

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

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

[Libpq-Fortran](https://github.com/ShinobuAmasaki/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つのクエリを実行しています。

```fortran
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`に記述する接続情報には、前のセクションで作成したロール・データベース・パスワードを使用します。

```fortran
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制約[^1]あり
- `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;`を問い合わせて、以下のような属性名の一覧を得られれば成功です。

```shell
$ 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]: この属性を空（`NULL`）にしてレコードを登録することができないという意味。

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

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

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

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

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

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

```fortran
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`の値に応じて切り替えて実行するようにしています。

##### ファイル出力

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

```fortran
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文と呼ばれるもので、テーブルにレコードを挿入します。

```fortran
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引数            | パラメータの個数だけの長さを持つゼロ値配列<br />（ここでのゼロはパラメータの型の解釈をサーバーに委ねるの意味） |
| 第5引数            | パラメータの個数だけの長さを持つ文字列型の配列               |

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


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

ここまでで必要な材料はすべて揃いました。メインプログラムを書いてシミュレーションを実行します。メインプログラムの全体は[GitHubのこちらのリンク](https://github.com/ShinobuAmasaki/ShinobuAmasaki.github.io/blob/main/sample-codes/sim-with-db/main.f90)から入手することができます。

```fortran
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で図を描画する方法について説明したいと思います[^2]。

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

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

以下が`plot.f90`の例です。

```fortran
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`に含まれるコマンドを実行しています。

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

```shell 
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スクリプトになります。

```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`）を実行すると次のような画像を出力することができます。

![plot9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2953366/b0c57942-2c94-caaf-d581-ee377d6602c9.png)


以上で、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