Mercurial > op > rk_sakamoti

annotate rungekutta.f90 @ 2:206e90e91191

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