※Canteraバージョン：3.0.0
Canteraは、Pythonの他にもC++やFORTRANなどの言語でも使用することができる。今回は、FortranプログラムでCanteraを使ってみる。
目次
ソースコードのコンパイル
CanteraをFortranで使うには、以前のC++の記事
と同様に、以下の公式サイトの方法で、ライブラリをインストールして使うか、
https://cantera.org/install/conda-install.html#sec-conda-development-interface
もしくは、以下のページのようにCanteraのソースコードをダウンロードして、ソースからコンパイルする。
https://cantera.org/install/compiling-install.html#sec-compiling
ただし、Canteraの既存のFortranインターフェースは、Canteraの一部の関数にしか対応しておらず、反応器Reactorの関数などは用意されていないようだ。
そこで、以下の例題
https://cantera.org/examples/fortran/index.html
にあるdemo_ftnlib.cppを参考に、C++でインターフェースとなる関数を作成して、これをFortranから呼び出す方法で試してみる。
今回も、CanteraのC++ライブラリは、WindowsのMinGW-W64のコンパイラーでソースからビルドしたものを用いた。FortranもMinGW-W64のgfortranでコンパイルした。
着火遅れ時間の計算
C++の時と同じく、上記の着火遅れ時間の計算プログラムをベースとして用いる。
Fortranプログラム
メインプログラムはFortranの方で、着火遅れを計算する関数ignitionDelayを呼び、着火遅れ時間を取得する。そして、計算結果として格納された時間と温度の時刻歴をCSVファイルに書き出す。
interface
double precision function ignitionDelay(mechFileName, temp, p)
character(len=*), intent(in) :: mechFileName
double precision, intent(in) :: temp, p
end function ignitionDelay
double precision function getStates(i, j)
integer, intent(in) :: i, j
end function getStates
subroutine getStatesSize(nr, nc)
integer, intent(out) :: nr, nc
end subroutine getStatesSize
end interface関数として、着火遅れを計算するignitionDelayと、計算結果（時間と温度）を取得するgetStatesを宣言する。また、サブルーチンとして結果配列のサイズを取得するgetStatesSizeを宣言する。
ignitionDelayは、反応メカニズムファイル名mechFileName、温度temp、圧力pを引数として受け取り、Canteraで計算して着火遅れ時間を返す。
getStatesは、時間と温度が格納された２次元配列のインデックスを引数として受け取り、時間または温度を返す。
getStatesSizeはサブルーチンで、結果の２次元配列のサイズが、nr、ncで取得できる。
これらの中身は、後述のC++のプログラムで実装されており、ここではinterface文で宣言している。
! input parameter
mechFileName = "LLNL_heptane_160.yaml"
temp = 1000.0
p = 1.3e6
! ignition delay C++ function
time_igd = ignitionDelay(trim(mechFileName), temp, p)
write(*, *) "Ignition Delay Time: ", time_igd * 1e6, " micro sec" ignitionDelay関数により、着火遅れ時間time_igdを取得する。
! output states
open(10, file="output.csv", status="replace")
write(10, "(a)") "time,temperature"
call getStatesSize(nr, nc)
do i = 0, nc-1
write(10, "(e0.8,a1,e0.8)") getStates(0, i), ",", getStates(1, i)
enddo
close(10)時間と温度をCSVファイルに出力する。サブルーチンgetStatesSizeを呼ぶと、nr、ncに２次元配列のサイズを取得できる。getStates(0,i)で時間、getStates(1,i)で温度を取得してファイルに書き出す。
C++プログラム
Fortranで呼び出す関数やサブルーチンの実装は、C++のプログラムで行う。もとになる内容は、上記のC++プログラミングの記事を参照。
static shared_ptr<Array2D> _states = nullptr;２次元配列Array2Dクラスの_statesに、計算結果の時間と温度を格納する。これは複数の関数で共有するため、冒頭で宣言しておく。ここでは、ポインタshared_ptrとする。
extern "C"C++の名前マングリング（name mangling）を避けるため、Fortranに渡す関数はextern "C"で記述する。
double ignitiondelay_(char* mechFileName, const double* temp, const double* p, ftnlen lenmech)着火遅れ時間を計算するignitionDelay関数を作成する。
関数名は小文字にして、最後に"_"をつけておく。
引数は、Fortranから参照で渡されるため、ポインタ引数とする。引数の最後のftnlen lenmechは、Fotran側では明示的に指定されていないが、文字列を引数にとる場合に追加される引数であり、文字列の長さが渡される。
std::string strMech = std::string(mechFileName, lenmech);引数で渡された文字列mechFileNameと文字列長さlenmechを使って、string型に変換している。
double getstates_(const int* i, const int* j)
{
if (_states == nullptr) {
return std::nan("");
}
size_t ist = *i;
size_t jst = *j;
if (ist >= 0 && ist < _states->nRows() && jst >= 0 && jst < _states->nColumns()) {
return _states->value(ist, jst);
} else {
return std::nan("");
}
}getStates関数は、インデックスを引数として受け取り、_statesに格納された時間または温度を返す。_statesが定義されていなかったり、インデックス範囲外であればNaNを返す。
void getstatessize_(int* nr, int* nc)
{
if (_states == nullptr) {
*nr = 0;
*nc = 0;
} else {
*nr = (int)_states->nRows();
*nc = (int)_states->nColumns();
}
}getStatesSizeは、２次元配列_statesのサイズを取得して、引数nr、ncの値を書き換える。
ここではあまりやっていないが、汎用的なライブラリとするには、きちんとエラー処理など行った方がよいだろう。
結果
Fortran、C++で各プログラムをコンパイル後、リンクして計算する。結果は、PythonやC++と同じになる。
自分のFortranプログラムに、Canteraの計算を組み込むことができる。