はじめに
ひょんなとこから Abaqus の UMAT をコーディングをすることになりました。言語は Fortran です。
Fortran なんて触ったこともないのでとりあえずググってみます。
すると Fortran の入門記事がいろいろと見つかるわけですが、どれもこれも最初に
implicit none を使え
と書いてあるわけです。先人がみな言うのですから、これはぜひとも従いたいところです。
ところが、Abaqus のユーザーサブルーチンのサンプルコードを見ても implicit none を使用しているものはほぼありません。いったいこれはどうしたことでしょう？
Abaqus のユーザーサブルーチンでは implicit none は使えないのでしょうか？
要点
結論から言うと、Abaqus のユーザーサブルーチンで implicit none は使えます。使えますが、ひと工夫必要です。
Abaqus のインターフェースやドキュメントは Fortran77 を引きずっていますが、Fortran90 以降のモダンな機能は使えます。
というわけで implicit none を使いつつ、module や function や type を使って読みやすいコードを書くとよいです。
方法
Abaqus のユーザーサブルーチンで implicit none を使うのが難しい理由は2つあります。
1つはユーザーサブルーチンの引数の型が暗黙に決まっていて分からないこと。
1つはインクルードファイル中で宣言されているパラメータで方が宣言されていないこと、です。
implicit none を宣言しない場合 Fortran では変数名で型が決まります。変数名が IJKLMN で始まっていれば INTEGER それ以外は REAL です。
が、ここにも落とし穴があります。REAL は単精度浮動小数点数を表す型ですが、Abaqus のような FEM の計算には単精度では足りないような……？
こうした謎を解くカギは、インクルードファイルにあります。ABA_PARAM.INC の中身を見てみましょう。これは Abaqus2018 の ABA_PARAM.INC です。
implicit real*8(a-h,o-z)
parameter (nprecd=2)
1行目に implicit real*8(a-h,o-z) とあります。デフォルトでは REAL 型となるところを REAL*8 型にせよということです。
REAL*8 は倍精度です。Abaqus は倍精度浮動小数点数で計算をしているのがこれでわかります。
というわけで、上記の1つ目の理由はこれで解決です。変数名が IJKLMN で始まっていれば INTEGER それ以外は REAL*8 ということになります。1
2行目の parameter (nprecd=2) ですが、これは定数の定義です。 n で始まっているのでこれは INTEGER です。
何に使っているのかわからないですが、おそらく Abaqus の組み込みサブルーチンなどで使うのでしょう。
implicit none を使うとこの定数の型が指定されていないためにエラーになります。これについては、インクルードファイルをあきらめて、コード中で都度定義することにします。
実際のコーディング
ここでは実例として弾性体の UMAT の例を見せます。
Fortran では自由記法と固定記法の二つの書き方がありますが、今回は固定記法を使います。理由はドキュメントやサンプルコードがほぼ固定記法になっているため、そこの変換の手間を避けるためです。
自由記法で不便がないなら自由記法を使う方がいいとは思います。
UMAT サブルーチン内に実装を書かかない
UMAT サブルーチン自体は Abaqus のドキュメントでインターフェースが提供されています。
これをなるべく触らないようにするため、 UMAT サブルーチン内では実装せず umat_impl サブルーチンを呼び出すだけとします。
umat_impl サブルーチンは関連する function やその他のサブルーチンとまとめるために implements モジュールの中に置くことにします。
そのように実装した UMAT が以下です。
UMAT は implicit none にしません。
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
use implements
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
4 JSTEP(4)
call umat_impl(STRESS,DDSDDE,DSTRAN,NTENS,PROPS,NPROPS)
RETURN
END
サブルーチン内でインクルードファイルの定数を定義する
呼び出すモジュール umat_impl はこのように定義します。
UMAT にわたってきた引数のうち、必要なものだけを引数として渡しています。弾性体だと必要なものは少ないです。
それらについては INTEGER だったり REAL8 だったり REAL8 の配列だったりするので UMAT のドキュメントやインターフェースと整合するように定義します。
module implements
implicit none
contains
subroutine umat_impl(STRESS,DDSDDE,DSTRAN,NTENS,PROPS,NPROPS)
! replace aba_param.inc for use implicit none
integer,parameter :: nprecd=2
integer :: NTENS,NPROPS
real*8 :: STRESS(NTENS),DDSDDE(NTENS,NTENS),DSTRAN(NTENS),
1 PROPS(NPROPS)
!
! Process
!
end subroutine umat_impl
end module implements
さらに integer,parameter :: nprecd=2 で ABA_PARAM.INC 内にあった定数を定義します。これについては umat_impl からさらにサブルーチンを呼ぶ場合はそのサブルーチン内でも定義する必要があるようです。
完成したコード
以上から、さらに内部の処理を umat_impl に書いたものがこれになります。type や function も使っていますが、これは使えるということを示す例として入れたもので正直使う必要はないです。
このコードは Tutorial: Write a simple UMAT in ABAQUS にあるコードをベースとして改変したものです。比較してみると何が変わったのかより分かりやすいです。
module implements
implicit none
type LAMBDA_struct
real*8 :: A,B,C
end type LAMBDA_struct
contains
function calc_LAMBDA(E,ANU)
type(LAMBDA_struct)::calc_LAMBDA
real*8::E,ANU
integer,parameter :: ONE=1.0D0, TWO=2.0D0
calc_LAMBDA%A=E/(ONE+ANU)/(ONE-TWO*ANU)
calc_LAMBDA%B=(ONE-ANU)
calc_LAMBDA%C=(ONE-TWO*ANU)
end function calc_LAMBDA
subroutine umat_impl(STRESS,DDSDDE,DSTRAN,NTENS,PROPS,NPROPS)
! replace aba_param.inc for use implicit none
integer,parameter :: nprecd=2
integer :: NTENS,NPROPS
real*8 :: STRESS(NTENS),DDSDDE(NTENS,NTENS),DSTRAN(NTENS),
1 PROPS(NPROPS)
integer :: I,J
type(LAMBDA_struct) :: LAMBDA
real*8 :: E,ANU
C ****************************************
C *** ELASTIC ELEMENT
C ****************************************
E=PROPS(1)
ANU=PROPS(2)
LAMBDA=calc_LAMBDA(E,ANU)
C ************************
DDSDDE=0.0D0
C ************************
DDSDDE(1,1)=(LAMBDA%A*LAMBDA%B)
DDSDDE(2,2)=(LAMBDA%A*LAMBDA%B)
DDSDDE(3,3)=(LAMBDA%A*LAMBDA%B)
DDSDDE(4,4)=(LAMBDA%A*LAMBDA%B)
DDSDDE(5,5)=(LAMBDA%A*LAMBDA%B)
DDSDDE(6,6)=(LAMBDA%A*LAMBDA%B)
DDSDDE(1,2)=(LAMBDA%A*ANU)
DDSDDE(1,3)=(LAMBDA%A*ANU)
DDSDDE(2,1)=(LAMBDA%A*ANU)
DDSDDE(2,3)=(LAMBDA%A*ANU)
DDSDDE(3,1)=(LAMBDA%A*ANU)
DDSDDE(3,2)=(LAMBDA%A*ANU)
C ************************
DO I=1,NTENS
DO J=1,NTENS
STRESS(I)=STRESS(I)+DDSDDE(I,J)*DSTRAN(J)
ENDDO
ENDDO
end subroutine umat_impl
end module implements
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
use implements
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
4 JSTEP(4)
call umat_impl(STRESS,DDSDDE,DSTRAN,NTENS,PROPS,NPROPS)
RETURN
END
最後に
Abaqus のユーザーサブルーチン界隈は、いまだに固定形式が当たり前だったりとプログラミングの手法がアップデートされていない印象があります。
固定形式はまだ我慢できますが、implicit none を使わないのは少しでもやりたい処理が複雑になると変数名がわかりにくくなる上に、変数名にかかわるバグで悩まされることになりいいことはありません。
Abaqus は implicit な変数定義を強要してくる非道なプロダクトですが、この記事を生かして implicit none を使いましょう。
ついでに module や function や subroutine や type を使ってコードにより意味を持たせ、処理を分割し、読みやすいコードを書くことで、ユーザーサブルーチンの実装を楽してもらえればと思います。
コードが複雑になればなるほど、コーディングの時間に占めるバグ取りの時間、つまりはコードを読み返す時間は増えます。多少書くのが面倒でも、読みやすさを優先してコーディングしたいものです。
参考文献
Tutorial: Write a simple UMAT in ABAQUS
编写ABAQUS子程序时启用implicit none
-
インクルードファイルは ABA_PARAM.INC 以外にもいくつかあります、Abaqus のバージョンによっても内容が変わるかもしれないのでそれぞれの環境で確認しましょう。 ↩