-
Notifications
You must be signed in to change notification settings - Fork 3
/
projector.f90
124 lines (90 loc) · 3.7 KB
/
projector.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
module projector
use numbers
use openmpi
use io
use parameters
use fieldio
use fftw
use vfield
complex(dpc), allocatable :: proj_bases(:,:,:,:,:), proj_origin(:, :, :, :)
real(dp), allocatable :: proj_results(:)
real(dp) :: proj_ekin
integer(i4) :: proj_ch, proj_uindex
character(1), parameter :: proj_uname = 'u'
character(255) :: proj_num_bases_str
character(255) :: proj_results_format
character(255) :: proj_dir
logical :: proj_written = .false.
contains
!==============================================================================
subroutine projector_init
integer :: i, un
write(proj_num_bases_str, *) num_proj_bases
proj_results_format = "(A2,"//i4_f//","//sp_f//","//dp_f//","//TRIM(proj_num_bases_str)//dp_f//")"
if (my_id == 0) allocate(proj_results(num_proj_bases))
if (subtract_origin) allocate(proj_origin(nx_perproc, ny_half, nz, 3))
allocate(proj_bases(nx_perproc, ny_half, nz, 3,num_proj_bases))
inquire(file='basis.in', exist=there)
if (.not. there) then
write(out,*) 'basis.in not found.'
flush(out)
stop
else
open(newunit=un,status='old',file='basis.in')
read(un,'(A)') proj_dir
close(un)
write(out, *) 'proj_dir:'
write(out, *) proj_dir
end if
! Read the average
if (subtract_origin) then
proj_uindex = 0
call projector_read(proj_origin)
end if
! Read the modes
do i=1, num_proj_bases
proj_uindex = i
call projector_read(proj_bases(:,:,:,:,i))
end do
! Create the projections file, replace it if it exists
if (my_id == 0) then
open(newunit=proj_ch,file='projections.gp',status='replace')
write(proj_ch, "(A2,"//i4_len//","//sp_len//","//dp_len//","//dp_len//")") "# ", "itime", "time", "proj_ekin", "projections"
proj_written = .true.
end if
end subroutine projector_init
!==============================================================================
subroutine projector_project(vfieldk)
complex(dpc), intent(in) :: vfieldk(:, :, :, :)
integer :: i
real(dp) :: inprod
complex(dpc) :: proj_vfieldk(nx_perproc, ny_half, nz, 3)
! Subtract the average
if (subtract_origin) then
proj_vfieldk(:,:,:,1:3) = vfieldk(:,:,:,1:3) - proj_origin(:,:,:,1:3)
else
proj_vfieldk(:,:,:,1:3) = vfieldk(:,:,:,1:3)
end if
do i = 1, num_proj_bases
call vfield_inprod(proj_vfieldk, proj_bases(:,:,:,1:3,i), inprod, .false.)
if (my_id == 0) proj_results(i) = inprod
end do
! Compute the energy
call vfield_norm2(proj_vfieldk, proj_ekin, .false.)
end subroutine projector_project
!==============================================================================
subroutine projector_write
! Write the projections
if (my_id == 0) then
write(proj_ch, TRIM(proj_results_format)) " ", itime, time, proj_ekin, proj_results
end if
end subroutine projector_write
!==============================================================================
subroutine projector_read(vfieldk)
complex(dpc), intent(out) :: vfieldk(:, :, :, :)
write(file_ext, "(i6.6)") proj_uindex
fname = TRIM(proj_dir)//'/'//proj_uname//'.'//file_ext
call fieldio_read(vfieldk)
end subroutine projector_read
!==============================================================================
end module projector