Mercurial > op > rk_sakamoti

comparison rungekutta.f90 @ 0:5bd0f2a15d2a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
import from https://qiita.com/sakamoti/items/de851e0d07aeef5be310
author "uncorrelated zombie" <uncorrelated@yahoo.co.jp>
date Tue, 23 Jul 2024 02:30:57 +0900
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:5bd0f2a15d2a
1 module rungekutta
2 use,intrinsic :: iso_fortran_env
3 private
4 !---
5 integer,parameter,public :: dpkind=real64
6 real(kind=dpkind),public,parameter :: PI=3.1415926535897932384626433_dpkind
7 !---
8 type,public :: rk
9 real(kind=dpkind) :: time=0e0 !<計算開始からの総経過時間
10 real(kind=dpkind) :: dt=1e-4 !<時間刻み
11 procedure(rhside),pass,pointer :: fun => null()
12 !<微分方程式の右辺を計算する関数
13 procedure(rksol),pass,pointer :: solver => null()
14 !<微分方程式の解法
15 contains
16 final :: final_rk !<構造体のファイナライズ処理
17 procedure :: setup !<数値積分の初期値、解法をセットする手続き
18 end type
19 !-------------
20 interface
21 !---
22 ! Right hand side function interface
23 function rhside(self,time,wk)
24 import :: dpkind,rk
25 class(rk),intent(inout) :: self !parameter
26 real(kind=dpkind) :: time
27 real(kind=dpkind),dimension(:) :: wk
28 real(kind=dpkind),dimension(size(wk)) :: rhside
29 end function
30 !---
31 ! Solver interface
32 subroutine rksol(self,tnext,wk)
33 import :: dpkind,rk
34 class(rk),intent(inout) :: self
35 real(kind=dpkind),intent(in) :: tnext !output time
36 real(kind=dpkind),dimension(:) :: wk
37 end subroutine
38 end interface
39 !============================--
40 contains
41 !---
42 pure elemental subroutine final_rk(self)
43 type(rk),intent(inout) :: self
44 if(associated(self%fun))nullify(self%fun)
45 if(associated(self%solver))nullify(self%solver)
46 end subroutine
47 !---
48 subroutine setup(self,fun,solver,dt)
49 class(rk),intent(inout) :: self
50 procedure(rhside) :: fun !微分方程式の右辺を計算する関数
51 character(len=*),intent(in) :: solver !ソルバーの選択
52 real(kind=dpkind),intent(in),optional :: dt !時間刻み幅
53 !時間刻みのセット
54 if(present(dt))then
55 self%dt=dt
56 endif
57 self%fun=>fun
58 !ソルバーの選択
59 select case(solver)
60 case("rk4")
61 !ルンゲクッタ法
62 self%solver=>rk4
63 case("rkf45")
64 !ルンゲクッタ・フェルベルグ法(刻み幅自動調整型)
65 !self%solver=>rkf45
66 case default
67 write(error_unit,*) "solver must one of the 'rk4, rkf45'"
68 endselect
69 end subroutine
70 !---
71 subroutine rk4(self,tn,wk)
72 class(rk),intent(inout) :: self
73 real(kind=dpkind),intent(in) :: tn !次に出力したい時間
74 real(kind=dpkind),dimension(:) :: wk
75 real(kind=dpkind),dimension(size(wk)) :: k1,k2,k3,k4
76 real(kind=dpkind) :: dt_tmp
77 real(kind=dpkind) :: dthalf
78 associate(time=>self%time,dt=>self%dt)
79 dt_tmp=dt
80 do
81 if(time+dt>tn)dt_tmp=tn-time
82 dthalf=dt_tmp/2e0
83 !!$OMP workshare
84 k1=self%fun(time,wk)
85 k2=self%fun(time+dthalf,wk+k1*dthalf)
86 k3=self%fun(time+dthalf,wk+k2*dthalf)
87 k4=self%fun(time+dt_tmp,wk+k3*dt_tmp)
88 !状態変数を更新
89 wk=wk+dt_tmp*(k1+2e0*k2+2e0*k3+k4)/6e0
90 !!$OMP end workshare
91 !1ステップ次の解が計算できたので、時間を進める
92 time=time+dt_tmp
93 !時間の計算範囲がtnを超えたら計算終了
94 if(time>=tn)then
95 exit
96 endif
97 enddo
98 end associate
99 end subroutine
100 end module