-
Notifications
You must be signed in to change notification settings - Fork 0
/
module_solver.f90
142 lines (99 loc) · 3.35 KB
/
module_solver.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
!
! module solver, to solver the equation of A and phi
!
module solver
use global
use auxiliary
use wfunction
implicit none
! public
complex(kind=k2),allocatable,save :: ka(:),kphi(:,:)
complex(kind=k2),allocatable,save :: suma0(:)
complex(kind=k2),allocatable,save :: sumphi0(:,:)
contains
!=========================runge-kutta 4 order ======================================
subroutine solver0(t,i)
implicit none
real(kind=k1),intent(in) :: t
integer, intent(in) :: i
integer(4) :: idicp,jdicp,itime
!
!==========================k1========================================================
!
call derivate(t,i)
! acoff
do idicp =1,n_string_ab
acoff(idicp) = acoff0(idicp) + ka(idicp)*dcmplx(0.5_k1)*prop%chstep
suma0(idicp) = acoff0(idicp) + ka(idicp)*dcmplx(1.0_k1/6.0_k1)*prop%chstep
enddo
! phi
do jdicp =1,system%nptot
do idicp =1, fedvr3d%nb_r*fedvr3d%nb_angle
phi(idicp,jdicp) = phi0(idicp,jdicp) + kphi(idicp,jdicp)*dcmplx(0.5_k1)*prop%chstep
sumphi0(idicp,jdicp) = phi0(idicp,jdicp) + kphi(idicp,jdicp)*dcmplx(1.0_k1/6.0_k1)*prop%chstep
enddo
enddo
!
!==========================k2========================================================
!
call derivate(t+0.5_k1*prop%hstep,i)
! acoff
do idicp =1,n_string_ab
acoff(idicp) = acoff0(idicp) + ka(idicp)*dcmplx(0.5_k1*prop%chstep)
suma0(idicp) = suma0(idicp) + ka(idicp)*dcmplx(prop%chstep/3.0_k1)
enddo
! phi
do idicp =1, fedvr3d%nb_r*fedvr3d%nb_angle
do jdicp =1,system%nptot
phi(idicp,jdicp) = phi0(idicp,jdicp) + kphi(idicp,jdicp)*dcmplx(0.5_k1*prop%chstep)
sumphi0(idicp,jdicp) = sumphi0(idicp,jdicp) + kphi(idicp,jdicp)*dcmplx(prop%chstep/3.0_k1)
enddo
enddo
!
!==========================k3========================================================
!
!!$ print*, ka
!!$ print*, acoff
!!$ print*, 'module_solver.f90'
!!$ stop
call derivate(t+0.5_k1*prop%hstep,i)
! acoff
do idicp =1,n_string_ab
acoff(idicp) = acoff0(idicp) + ka(idicp)*dcmplx(prop%chstep)
suma0(idicp) = suma0(idicp) + ka(idicp)*dcmplx(prop%chstep/3.0_k1)
enddo
! phi
do idicp =1, fedvr3d%nb_r*fedvr3d%nb_angle
do jdicp =1,system%nptot
phi(idicp,jdicp) = phi0(idicp,jdicp) + kphi(idicp,jdicp)*dcmplx(prop%chstep)
sumphi0(idicp,jdicp) = sumphi0(idicp,jdicp) + kphi(idicp,jdicp)*dcmplx(prop%chstep/3.0_k1)
enddo
enddo
!
!==========================k4========================================================
!
call derivate(t+prop%hstep,i)
! acoff
do idicp =1,n_string_ab
acoff(idicp) = suma0(idicp) + ka(idicp)*dcmplx(prop%chstep/6.0_k1)
enddo
! phi
do idicp =1, fedvr3d%nb_r*fedvr3d%nb_angle
do jdicp =1,system%nptot
phi(idicp,jdicp) = sumphi0(idicp,jdicp) + kphi(idicp,jdicp)*dcmplx(prop%chstep/6.0_k1)
enddo
enddo
return
end subroutine solver0
end module solver
subroutine configure_solver
use global
use auxiliary
use solver
implicit none
allocate(ka(n_string_ab))
allocate(kphi(fedvr3d%nb_r*fedvr3d%nb_angle,system%nptot))
allocate(suma0(n_string_ab))
allocate(sumphi0(fedvr3d%nb_r*fedvr3d%nb_angle,system%nptot))
return
end subroutine configure_solver