-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathesvm.f90
177 lines (177 loc) · 6.76 KB
/
esvm.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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! ElectroStatic Vlasov-Maxwell (ESVM) code !!
!! !!
!! Copyright © 2015 Michaël J TOUATI !!
!! !!
!! This file is part of ESVM. !!
!! !!
!! ESVM is free software: you can redistribute it and/or modify !!
!! it under the terms of the GNU General Public License as published !!
!! by the Free Software Foundation, either version 3 of the License, !!
!! or (at your option) any later version. !!
!! !!
!! ESVM is distributed in the hope that it will be useful, !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!
!! GNU General Public License for more details. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with ESVM. If not, see <https://www.gnu.org/licenses/>. !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Initial commit written by Michaël J TOUATI - Dec. 2015
program ESVM
!
use acuracy
use constants
use input
use diagnostics
use library
!
implicit none
real(PR) :: time, d_t
real(PR), dimension(:), allocatable :: x, vx
real(PR), dimension(:,:), allocatable :: f_n, f_np1
real(PR), dimension(:), allocatable :: n_e, j_e, v_e, vT_e
real(PR), dimension(:), allocatable :: E_x_n, E_x_np1, phi_n
integer :: N_t, i, l
real(PR) :: f_max, flux_im1, flux_ip1
real(PR), dimension(:), allocatable :: dU_K, dU_T, dU_E
real(PR) :: U_K, U_T, U_E
real(PR) :: vm1, am1
logical :: test_positivity
logical :: save_energies, save_results
real(PR) :: timer_start, timer_finish
real(PR) :: t_drive, dt_diag_nrj
real(PR) :: CPUtime, clock
integer :: clock_start, clock_finish
integer :: rate
!
call system_clock(count_rate=rate)
call system_clock(clock_start)
call cpu_time(timer_start)
!
call read_init_parameters()
!
call system('mkdir -p results/'//trim(simu))
!
allocate(x(-1:N_x+2),vx(-1:N_vx+2))
allocate(f_n(-1:N_x+2,-1:N_vx+2),f_np1(-1:N_x+2,-1:N_vx+2))
allocate(n_e(-1:N_x+2),j_e(-1:N_x+2),v_e(-1:N_x+2),vT_e(-1:N_x+2))
allocate(phi_n(-1:N_x+2),E_x_n(-1:N_x+2),E_x_np1(-1:N_x+2))
allocate(dU_K(1:N_x),dU_T(1:N_x),dU_E(1:N_x))
!
call GRID(N_x, N_vx, d_x, d_vx, &
& x_min, vx_min, x(-1:N_x+2), vx(-1:N_vx+2))
!
call INIT_VAR(N_x, N_vx, f_n, f_np1, &
& n_e, j_e, v_e, vT_e, &
& E_x_np1, E_x_n , phi_n, &
& dU_K, dU_T, dU_E, &
& U_K, U_T, U_E, time, N_t, &
& test_positivity, save_results)
!
call INIT_SIMU(b_cond, perturb, A_pert, k_pert, v_d, &
& N_x, N_vx, x(-1:N_x+2), vx(-1:N_vx+2), f_n)
!
call INIT_DIAG()
!
t_drive = 6._PR * pi / omega_pert
dt_diag_nrj = 0.1_PR
if (perturb == 2) dt_diag_nrj = dt_diag_nrj / omega_pert
!
do while (time.lt.L_t)
!
f_max = maxval(maxval(f_n(:,:),dim=2))
!
call DENSITIES(N_x, N_vx, d_vx, &
& vx(-1:N_vx+2), f_n(-1:N_x+2,-1:N_vx+2),&
& n_e(-1:N_x+2), j_e(-1:N_x+2),&
& v_e(-1:N_x+2), vT_e(-1:N_x+2))
! Landau damping test-case
if (perturb == 2) then
if ( time < t_drive ) then
call DRIVE(N_x, d_t, time, x, &
& A_pert, omega_pert, k_pert, &
& E_x_n, E_x_np1, phi_n)
else
call MAXWELL_SOLVER(maxwell, b_cond, &
& N_x, N_t, d_t, d_x, &
& j_e, n_e, E_x_n, E_x_np1, phi_n)
end if
! All other cases
else
call MAXWELL_SOLVER(maxwell, b_cond, &
& N_x, N_t, d_t, d_x, &
& j_e, n_e, E_x_n, E_x_np1, phi_n)
end if
!
call ENERGIES(N_x, d_x, &
& n_e, v_e, vT_e, E_x_n, &
& dU_K, dU_T, dU_E, &
& U_K, U_T, U_E)
!
d_t = cfl*(0.5_PR/((maxval(vx(1:N_vx))/d_x)+(maxval(abs(E_x_n(1:N_x)))/d_vx)))
vm1 = d_t / d_x
am1 = d_t / d_vx
!
!$omp PARALLEL DO DEFAULT(NONE) &
!$omp& SHARED(N_vx, N_x, d_t, d_x, d_vx, f_n, f_np1, vx, E_x_n) &
!$omp& SHARED(vm1, am1, test_positivity, scheme, b_VL, f_max) &
!$omp& PRIVATE(l,i,flux_im1,flux_ip1) COLLAPSE(2)
do l=1,N_vx,1
do i = 1,N_x,1
call fluxes(scheme, b_VL, vx(l), f_max, d_t, d_x, &
& f_n(i-2,l), f_n(i-1,l), f_n(i,l), f_n(i+1,l), f_n(i+2,l), &
& flux_im1, flux_ip1)
f_np1(i,l) = f_n(i,l) - ( vm1 * (flux_ip1 - flux_im1) )
!
call fluxes(scheme, b_VL, -E_x_n(i), f_max, d_t, d_vx,&
& f_n(i,l-2), f_n(i,l-1), f_n(i,l), f_n(i,l+1), f_n(i,l+2), &
& flux_im1, flux_ip1)
f_np1(i,l) = f_np1(i,l) - ( am1 * (flux_ip1 - flux_im1) )
if ((f_np1(i,l).lt.0._PR).and.(test_positivity.eqv..false.)) test_positivity = .true.
end do
end do
!$omp END PARALLEL DO
!
call FE_BOUNDARIES(b_cond, N_x, N_vx, f_np1(-1:N_x+2,-1:N_vx+2))
!
save_energies = (mod(time,dt_diag_nrj).lt.d_t).and.(time.ge.d_t)
save_energies = save_energies.or.(N_t.eq.1)
save_energies = save_energies.or.((L_t-time).le.d_t)
!
if (save_energies.eqv..true.) then
call DIAG_ENERGY(time, U_K, U_T, U_E)
end if
!
save_results = (mod(time,dt_diag).lt.d_t).and.(time.ge.d_t)
save_results = save_results.or.(N_t.eq.1)
save_results = save_results.or.((L_t-time).le.d_t)
!
if (save_results.eqv..true.) then
call DIAG(N_t, time, N_x, x, N_vx, vx, &
& test_positivity, U_K, U_T, U_E, &
& f_n, n_e, E_x_n, j_e, v_e, vT_e, phi_n)
end if
!
call INIT_NEXT_STEP(N_x, N_vx, f_n, f_np1, &
& E_x_n, E_x_np1, &
& N_t, d_t, time)
end do
!
deallocate(x, vx, f_n, f_np1)
deallocate(n_e, j_e, v_e, vT_e)
deallocate(E_x_n, E_x_np1, phi_n)
deallocate(dU_K, dU_T, dU_E)
!
call system_clock(clock_finish)
call cpu_time(timer_finish)
clock = real(clock_finish - clock_start,PR)/real(rate,PR)
clock = clock/3.6e3_PR
CPUtime = (timer_finish - timer_start)/3.6e3_PR
!
call CLOSE_DIAG(CPUtime, clock)
!
end program ESVM