能書き
以前、二次元ラプラス方程式をキャラクタ・グラフィックスで表示してみましたが、それを内容そのままで出力形式をアニメーション gif にしてみます。
アニメーション gif のルーチンは、 fortran66 さんから借りてくることにします。初期化と終了処理を除けば、以前のキャラクター出力の所を少しいじるだけで直せます。色の設定が難しいですね。
ラプラス方程式の解は、微分方程式を反復法で求めています。ここで mask 付きの forall 構文を用いることで、微分方程式の反復解法の記述と境界条件の記述とを分離することができます。
プログラムでは往々にして、解法の記述と境界条件の記述とが渾然一体となってしまうので、この方式だとスッキリ気分を味わえます。
出力結果
外側のパイプが 100 度、内側のパイプが 0 度、パイプ間の媒質が 20 度という初期値で計算しています。図では温度を 10 度刻みで色の濃淡で表わしています。
同心円の場合
! define shape : concentric circle
mask = .false.
forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mask(i, j) = .true.
forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 50**2) mask(i, j) = .false.
! initial value
mesh = 100.0
forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mesh(i, j) = 20.0
forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 50**2) mesh(i, j) = 0.0
中心のずれた円の場合
! define shape : eccentric pipe
mask = .false.
forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mask(i, j) = .true.
forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 351)**2 < 50**2) mask(i, j) = .false.
! initial value
mesh = 100.0
forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mesh(i, j) = 20.0
forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 351)**2 < 50**2) mesh(i, j) = 0.0
内側のパイプが矩形の場合
! define shape : circle and square
mask = .false.
forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2 ) mask(i, j) = .true.
forall (i = 1:n, j = 1:n, abs(i - 351) < 75 .and. abs(j - 351) < 75) mask(i, j) = .false.
! initial value
mesh = 100.0
forall (i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2 ) mesh(i, j) = 20.0
forall (i = 1:n, j = 1:n, abs(i - 351) < 75 .and. abs(j - 351) < 75) mesh(i, j) = 0.0
プログラム・リスト
プログラムでは、Fortran2008 の機能を使っています。forall(integer:: ) の構文で局所変数を使っているので、Intel Fortran の Ver.19 beta 版でないと動きません。（integer を普通に変数宣言すれば F95 形式になります。)なお、コンパイラのオプションで、semantics を Fortran2003 対応にしておかないと、微妙なメモリー・エラーが出たり出なかったりします。
Fortran2008 は block 構文が導入されて、局所変数が使えるだけでなく、use 文も使えるようになりました。サブルーチン化するほどでもないが局所的にまとめておきたい、たとえば今回の様なグラフィックスのような目的に便利な機能です。
!
! animation gif modules borrowed from http://fortran66.hatenablog.com/entry/2018/05/09/014127
!
module m_gif_types
use, intrinsic :: iso_fortran_env
implicit none
private
public :: t_gif_header, t_image_block, t_graphic_control_extension, t_application_extension
! constructors (interface)
interface t_gif_header
procedure :: init_gif_header
end interface
interface t_image_block
procedure :: init_image_block
end interface
interface t_graphic_control_extension
procedure :: init_graphic_control_extension
end interface
interface t_application_extension
procedure :: init_application_extension
end interface
!
type :: t_gif_header
sequence
character(3) :: signature = 'GIF'
character(3) :: version = '89a' ! '87a'
integer(int16) :: width
integer(int16) :: height
integer(int8) :: pck = int(B'10001010', int8) !1:Global Col. 3:Color Res. 1:Sort Flg. 3:Size of Global Col.
integer(int8) :: background_color_index = 0
integer(int8) :: pixel_aspect_ratio = 0
end type
type :: t_image_block
sequence
integer(int8) :: image_separator = int(Z'2C', int8)
integer(int16) :: image_left_position = 00
integer(int16) :: image_top_position = 00
integer(int16) :: image_width
integer(int16) :: image_height
integer(int8) :: pck = int(B'00000000', int8)
! local color table
!integer(int8) :: LZW_minimum_code_size
!integer(int8) :: block_size
!integer(int8) :: image_data (:)
!integer(int8) :: block_terminator = int(Z'00', int8)
end type t_image_block
type :: t_graphic_control_extension
sequence
integer(int8) :: extention_introducer = int(Z'21', int8)
integer(int8) :: graphic_control_label = int(Z'F9', int8)
integer(int8) :: lock_size = int(Z'04', int8)
integer(int8) :: pck = int(Z'00', int8)
integer(int16) :: delay_time = 50 ! msec
integer(int8) :: transparent_color_index = 0
integer(int8) :: block_terminator = int(Z'00', int8)
end type t_graphic_control_extension
type :: t_application_extension ! for animation gif
sequence
integer(int8) :: extension_intrducer = int(Z'21', int8)
integer(int8) :: extension_label = int(Z'FF', int8)
integer(int8) :: block_size_01 = int(Z'0B', int8)
character(len = 8) :: application_identifier = 'NETSCAPE'
character(len = 3) :: application_authentication_code = '2.0'
integer(int8) :: block_size_02 = int(Z'03', int8) ! 0:block terminator
integer(int8) :: n = int(Z'01', int8) !application_data
integer(int16) :: nloop = 0 ! 0:unlimited loop
integer(int8) :: block_terminator = int(Z'00', int8)
end type t_application_extension
contains
! constructors (implementation)
pure type(t_gif_header) function init_gif_header(nx, ny, ngcol) result(res)
integer, intent(in) :: nx, ny, ngcol
res%width = nx
res%height = ny
res%pck = iand(res%pck, int(B'11111000', int8)) + int(mod(ngcol, 8), int8) !least 3bits:size of global color:2**(ngcol + 1) 2..256
end function init_gif_header
pure type(t_image_block) function init_image_block(nx, ny) result(res)
integer, intent(in) :: nx, ny
res%image_width = nx
res%image_height = ny
end function init_image_block
pure type(t_graphic_control_extension) function init_graphic_control_extension(it) result(res)
integer, intent(in) :: it
res%delay_time = int(it, int16) ! msec
end function init_graphic_control_extension
pure type(t_application_extension) function init_application_extension(nloop) result(res)
integer, intent(in) :: nloop
res%nloop = int(nloop, int16)
end function init_application_extension
end module m_gif_types
module m_gif_lzw
use, intrinsic :: iso_fortran_env
use :: m_gif_types
implicit none
private
public :: encoder
type :: t_enc
integer :: kbits = 0, id = 0
end type
contains
subroutine enc(irgb, nbits, code) ! gif LZW
integer, intent(in) :: irgb(:), nbits
type (t_enc), intent(out) :: code(:)
integer :: dict(0:2**12 - 1) ! 0:4095
integer :: kbits, kd, m0 ! <== m0 state variable
integer :: ip, ic, ienc
call clear_dict() ! clear dictionary
ic = 1
code(ic) = t_enc(kbits, 2**nbits) ! clear_code
m0 = irgb(1)
do ip = 2, size(irgb)
call subenc(ienc, irgb(ip))
if (ienc /= -1) then ! new word
ic = ic + 1
code(ic) = t_enc(kbits, ienc)
end if
if (kd == 2**12 - 1) then ! dictionary full
ic = ic + 1
code(ic) = t_enc(kbits, 2**nbits) ! clear_code
call clear_dict() ! clear dictionary
end if
end do
ic = ic + 1
code(ic) = t_enc(kbits, m0)
ic = ic + 1
code(ic) = t_enc(kbits, dict(2**nbits + 1)) ! end_code
return
contains
subroutine subenc(ienc, m1)
integer, intent(out) :: ienc
integer, intent(in ) :: m1
integer :: k, id
k = ishft(m0 + 1, 16) + m1 ! m0 + 1 to avoid 00...0 degeneracy
id = findloc(dict, k, dim = 1) ! dictionary is 0-based
if (id == 0) then ! not found in the dictionary
kd = kd + 1
dict(kd) = k
if (kd == 2**kbits + 1) kbits = kbits + 1
ienc = m0
m0 = m1
else ! found in the dictionary
ienc = -1
m0 = id - 1 ! because dictionary is 0-based, m0 must be shifted by 1
end if
end subroutine subenc
subroutine clear_dict()
kbits = nbits + 1
kd = 2**nbits + 1 ! dictinary defined
dict = 0
forall(integer :: i = 0:kd) dict(i) = i ! dict(nbit):clear code; dict(nbit+1):end code
end subroutine clear_dict
end subroutine enc
subroutine encoder(irgb, nbits, icode)
integer, intent(in ) :: irgb(:), nbits
integer(int8), allocatable, intent(out) :: icode(:)
type(t_enc), allocatable :: code(:)
integer :: i, j, k, ib, ic, nb
allocate(code(size(irgb) + 100)) ! ? large enough ?
call enc(irgb, nbits, code) ! compress color code to LZW code
nb = ceiling(sum(code(:)%kbits) / 8.0) ! required bytes
allocate(icode(nb))
icode = 0 ! pack LZW code to bit string
k = 0
do i = 1, size(code)
do j = 1, code(i)%kbits
if (btest(code(i)%id, j - 1)) then
ic = k / 8 + 1
ib = mod(k, 8)
icode(ic) = ibset(icode(ic), ib)
end if
k = k + 1
end do
end do
end subroutine encoder
end module m_gif_lzw
module m_anime_gif
use, intrinsic :: iso_fortran_env
use :: m_gif_types
use :: m_gif_lzw
implicit none
private
public :: t_gif
interface t_gif
procedure :: init_gif
end interface
type :: t_gif
integer :: iw, nbits
integer, allocatable :: irgb2(:, :)
contains
procedure :: wr_gif_img, gif_close
end type t_gif
contains
impure type(t_gif) function init_gif(fn, nx, ny, n_global_color, global_color_table, loop) result(this)
character(*), intent(in) :: fn
integer, intent(in) :: nx, ny, n_global_color, global_color_table(:)
integer, intent(in), optional :: loop
type(t_gif_header) :: head
type(t_application_extension) :: anime
integer :: nloop = 0
if (present(loop)) nloop = loop
associate(iw => this%iw, nbits => this%nbits)
nbits = max(2, n_global_color + 1) ! 1..8 -> 2,2,3,4,5,6,7,8
allocate(this%irgb2(nx, ny))
head = t_gif_header(nx, ny, n_global_color)
anime = t_application_extension(nloop)
open(newunit = iw, file = fn, access = 'stream')
write(iw) head
write(iw) int(global_color_table, int8)
write(iw) anime
end associate
end function init_gif
subroutine wr_gif_img(this, itime)
class(t_gif), target, intent(in) :: this
integer, intent(in) :: itime
type(t_graphic_control_extension), allocatable :: ext
type(t_image_block), allocatable :: img
integer(int8), allocatable :: icode(:)
integer :: i, nbyte
integer, pointer :: irgb(:)
associate(iw => this%iw, nbits => this%nbits, irgb2 => this%irgb2)
ext = t_graphic_control_extension(itime)
img = t_image_block(size(irgb2, 1), size(irgb2, 2))
irgb(1:size(irgb2)) => irgb2 ! 1D => 2D
call encoder(irgb, nbits, icode) ! color code to GIF LZW code
write(iw) ext ! write bit strings to file
write(iw) img
write(iw) int(nbits, int8) ! LZW_minimum_code_size = 2,2,3,4,5,6,7,8bits
do i = 1, size(icode), 255
nbyte = min(i + 255, size(icode)) - i
write(iw) achar(nbyte)
write(iw) icode(i:i+nbyte-1)
end do
write(iw) achar(00) ! block_terminator
end associate
end subroutine wr_gif_img
subroutine gif_close(this)
class(t_gif), intent(in out) :: this
associate(iw => this%iw)
write(iw) achar(int(Z'3B'))
close(iw)
end associate
end subroutine gif_close
end module m_anime_gif
program Heat
implicit none
integer, parameter :: n = 501
real :: mesh(0:n + 1, 0:n + 1)
logical :: mask(n, n)
integer :: iter
! define shape : eccentric pipe
mask = .false.
forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mask(i, j) = .true.
forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 351)**2 < 50**2) mask(i, j) = .false.
! initial value
mesh = 100.0
forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 251)**2 < 250**2) mesh(i, j) = 20.0
forall (integer:: i = 1:n, j = 1:n, (i - 251)**2 + (j - 351)**2 < 50**2) mesh(i, j) = 0.0
! solve Laplace's equation by iterative method
block
use m_anime_gif
type(t_gif), allocatable :: fig
integer :: n_global_color ! 0..7 ! 2..256 colors (2^(n_global_color + 1))
integer, allocatable :: global_color_table(:)
n_global_color = 7
allocate(global_color_table(3 * 2**(n_global_color + 1)))
global_color_table = 0
forall (integer :: i = 0:254) global_color_table(3*i+1:3*i+3) = mod(2 * i, 255)
fig = t_gif('heat.gif', n+2, n+2, n_global_color, global_color_table)
!
do iter = 0, 100000
if (mod(iter, 5000) == 0) then ! gif animation
fig%irgb2 = 10 * int(mesh / 10.0)
call fig%wr_gif_img(20)
end if
forall (integer:: i = 1:n, j = 1:n, mask(i, j)) &
mesh(i, j) = (mesh(i - 1, j) + mesh(i + 1, j) + mesh(i, j - 1) + mesh(i, j + 1)) / 4
end do
call fig%gif_close()
deallocate(fig)
end block
end program Heat