Newer
Older
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
143
144
145
146
147
148
149
150
SUBROUTINE moments_eq_rhs
USE basic
USE time_integration
USE array
USE fields
USE space_grid
USE model
use prec_const
IMPLICIT NONE
INTEGER:: ip, iz
REAL(dp) :: ip_dp
REAL(dp) :: tmp, moml, sqrtT
do ip = ips, ipe
ip_dp = real(ip,dp) ! From int to double (compute once)
do iz = izs,ize ! Compute coefficients
sqrtT = sqrt_exp_temp(iz)
! N_e^{p} term
moml = moments(ip,iz,updatetlevel)
moments_Apl(iz) = -nu*ip_dp*moml ! Collisional damping from Lenard-Bernstein Operator
moments_Bpl(iz) = -moml
moments_Cpl(iz) = -sqrtT*vpar_n(iz)*INVSQRT2*moml
moments_Dpl(iz) = 0._dp
moments_Epl(iz) = -ip_dp*0.5_dp*moml
moments_Fpl(iz) = -vpar_n(iz)*2._dp*SQRT2*ip_dp*moml
! moments_Fpl(iz) = -sqrtT*vpar_n(iz)*SQRT2*ip_dp*moml
moments_Gpl(iz) = 0._dp
moments_Hpl(iz) = -sqrtT*SQRT2*(ip_dp+1._dp)*moml
! N_e^{p+1} term
if (ip+1 .le. ipe) then
moml = moments(ip+1,iz,updatetlevel)
moments_Cpl(iz) = moments_Cpl(iz) -sqrtT*sqrt(ip_dp+1._dp)*0.5_dp*moml
moments_Fpl(iz) = moments_Fpl(iz) -ip_dp*sqrt(ip_dp+1._dp)*moml
! moments_Fpl(iz) = moments_Fpl(iz) -sqrtT*ip_dp*sqrt(ip_dp+1._dp)*0.5_dp*moml
endif
! N_e^{p-1} term
if (ip-1 .ge. ips) then
moml = moments(ip-1,iz,updatetlevel)
if (ip .ne. 3) then
moments_Apl(iz) = moments_Apl(iz) -nu*SQRT2*sqrt(ip_dp)*moml*vpar_n(iz) ! Collisions
endif
moments_Cpl(iz) = moments_Cpl(iz) -sqrtT*sqrt(ip_dp)*0.5_dp*moml
moments_Dpl(iz) = moments_Dpl(iz) +sqrt(ip_dp)/sqrtT*moml
moments_Epl(iz) = moments_Epl(iz) -vpar_n(iz)*sqrt(ip_dp*0.5_dp)*moml
moments_Fpl(iz) = moments_Fpl(iz) -sqrt(ip_dp)*(2._dp*ip_dp-1._dp+2._dp*vpar_n(iz)*vpar_n(iz))*moml
! moments_Fpl(iz) = moments_Fpl(iz) -sqrtT*sqrt(ip_dp)*0.5_dp*(2._dp*ip_dp-1._dp+2._dp*vpar_n(iz)*vpar_n(iz))*moml
moments_Gpl(iz) = moments_Gpl(iz) -SQRT2*sqrt(ip_dp)*moml
moments_Hpl(iz) = moments_Hpl(iz) -sqrtT*vpar_n(iz)*2._dp*sqrt(ip_dp)*moml
endif
! N_e^{p-2} term
if (ip-2 .ge. ips) then
moml = moments(ip-2,iz,updatetlevel)
moments_Epl(iz) = moments_Epl(iz) -sqrt(ip_dp*(ip_dp-1._dp))*0.5_dp*moml
moments_Fpl(iz) = moments_Fpl(iz) -2._dp*sqrt(2._dp*ip_dp*(ip_dp-1._dp))*moml
! moments_Fpl(iz) = moments_Fpl(iz) -sqrtT*sqrt(2._dp*ip_dp*(ip_dp-1._dp))*moml
moments_Hpl(iz) = moments_Hpl(iz) -sqrtT*sqrt(2._dp*ip_dp*(ip_dp-1._dp))*moml
endif
! N_e^{p-3} term
if (ip-3 .ge. ips) then
moments_Fpl(iz) = moments_Fpl(iz) -sqrt(ip_dp*(ip_dp-1._dp)*(ip_dp-2._dp))*moments(ip-3,iz,updatetlevel)
! moments_Fpl(iz) = moments_Fpl(iz) -sqrtT*sqrt(ip_dp*(ip_dp-1._dp)*(ip_dp-2._dp))*0.5_dp*moments(ip-3,iz,updatetlevel)
else if (ip .eq. ips) then ! ip-3=0 so N_e^0 is equal to 1
moments_Fpl(iz) = moments_Fpl(iz) -SQRT3*SQRT2
! moments_Fpl(iz) = moments_Fpl(iz) -sqrtT*SQRT3*SQRT2*0.5_dp
endif
end do
select case (gradient_scheme)
! case('up2') ! Upwind scheme, check the sign of coefficients
! do iz = izs,ize
! if (moments_Cpl(iz)>0) then
! tmp = moments_Cpl(iz)*thetaz(iz)
! else
! tmp = moments_Cpl(iz)*thetazm(iz)
! endif
! if (moments_Dpl(iz)>0) then
! tmp = tmp + moments_Dpl(iz)*phiz(iz)
! else
! tmp = tmp + moments_Dpl(iz)*phizm(iz)
! endif
! if (moments_Fpl(iz)>0) then
! tmp = tmp + moments_Fpl(iz)*tempz(iz)
! else
! tmp = tmp + moments_Fpl(iz)*tempzm(iz)
! endif
! if (moments_Hpl(iz)>0) then
! tmp = tmp + moments_Hpl(iz)*vparz(iz)
! else
! tmp = tmp + moments_Hpl(iz)*vparzm(iz)
! endif
! if (vpar(iz,updatetlevel)<0) then ! I_p^l
! moments_Ipl = -sqrtT*SQRT2*vpar(iz,updatetlevel)*momentsz(ip,iz)
! else
! moments_Ipl = -sqrtT*SQRT2*vpar(iz,updatetlevel)*momentszm(ip,iz)
! endif
! if (ip+1 .le. ipe) moments_Ipl = moments_Ipl -sqrtT*sqrt(ip_dp+1._dp)*momentszm(ip+1,iz) ! coefficient always negative
! if (ip-1 .ge. ips) moments_Ipl = moments_Ipl -sqrtT*sqrt(ip_dp)*momentszm(ip-1,iz) ! coefficient always negative
! moments_rhs(ip,iz,updatetlevel) = tmp &
! +moments_Apl(iz) &
! +moments_Bpl(iz)*theta_rhs(iz,updatetlevel) &
! +moments_Epl(iz)*temp_rhs(iz,updatetlevel) &
! +moments_Gpl(iz)*vpar_rhs(iz,updatetlevel) &
! +moments_Ipl
! end do
case('fa2','fd4','we4')
do iz = izs,ize
moments_Ipl = -sqrtT*SQRT2*vpar_n(iz)*momentsz(ip,iz)
if (ip+1 .le. ipe) moments_Ipl = moments_Ipl -sqrtT*sqrt(ip_dp+1._dp)*momentsz(ip+1,iz) ! coefficient always negative
if (ip-1 .ge. ips) moments_Ipl = moments_Ipl -sqrtT*sqrt(ip_dp)*momentsz(ip-1,iz) ! coefficient always negative
moments_rhs(ip,iz,updatetlevel) = diff_moments*momentszz(ip,iz) & ! numerical diffusion added for stability
+moments_Apl(iz) &
+moments_Bpl(iz)*theta_rhs(iz,updatetlevel) &
+moments_Cpl(iz)*thetaz(iz) &
+moments_Dpl(iz)*phiz(iz) &
+moments_Epl(iz)*temp_rhs(iz,updatetlevel) &
+moments_Fpl(iz)*sqrt_exp_tempz(iz) &
+moments_Gpl(iz)*vpar_rhs_n(iz) &
+moments_Hpl(iz)*vparz_n(iz) &
+moments_Ipl
! moments_rhs(ip,iz,updatetlevel) = 0._dp ! debug
end do
end select
end do
END SUBROUTINE moments_eq_rhs