-
Notifications
You must be signed in to change notification settings - Fork 0
/
init.f90
216 lines (196 loc) · 6.35 KB
/
init.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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
module init
use workspace
use io
contains
subroutine initCond
integer :: f, n
integer, dimension(:), allocatable :: seed
if (RANK .eq. 0) then
write (6, "(a)") 'Applying initial condition...'
end if
TIME = 0.0d0
if (RANK .eq. 0) then
write (6, "(a)") 'Seeding RNG'
end if
call RANDOM_SEED(size=n)
allocate (seed(n))
seed = RANK + RSEED
call RANDOM_SEED(PUT=seed)
if (initialCondType .eq. 1) then
do f = 1, FLUIDS
call makeRandomPhase(WS%FLUID(f))
end do
else if (initialCondType .eq. 2) then
call makeNonEquibPhase
else if (initialCondType .eq. 3) then
call read_wf_file(ICRfilename)
else if (initialCondType .eq. 4) then
do f = 1, FLUIDS
call makeRandomDens(WS%FLUID(f), 0.05d0)
end do
else if (initialCondType .eq. 5) then
include 'ic.in'
else
do f = 1, FLUIDS
call makeConst(WS%FLUID(f), (1.0d0, 0.0d0))
end do
end if
end subroutine
subroutine makeConst(field, c)
implicit none
class(computational_field) :: field
integer :: i, j, k
complex*16 :: c
!$OMP parallel do private (i,j,k) collapse(3)
do k = sz, ez
do j = sy, ey
do i = sx, ex
field%GRID(i, j, k) = c
end do
end do
end do
!$OMP end parallel do
end subroutine
subroutine makeRandomPhase(field)
implicit none
type(fluid_field) :: field
double precision :: r1
integer :: i, j, k
if (RANK .eq. 0 .and. ANY(OMEGAALL .ne. 0.0d0)) then
write (6, '(a,i3)') "Imposing a random phase IC on FIELD ", field%field_number
end if
!$OMP parallel do private (i,j,k) collapse(3)
do k = sz, ez
do j = sy, ey
do i = sx, ex
call random_number(r1)
field%GRID(i, j, k) = 1.d0*exp(2.0*PI*EYE*r1)
end do
end do
end do
!$OMP end parallel do
end subroutine
subroutine makeRandomDens(field, s)
implicit none
type(fluid_field) :: field
double precision :: r1, s
integer :: i, j, k
if (RANK .eq. 0) then
write (6, *) ' Imposing a random density perturbation on FIELD ', field%field_number
end if
!$OMP parallel do private (i,j,k) collapse(3)
do k = sz, ez
do j = sy, ey
do i = sx, ex
call random_number(r1)
field%GRID(i, j, k) = 1.0d0 + r1*s - s/2.0d0
end do
end do
end do
!$OMP end parallel do
end subroutine
subroutine makeNonEquibPhase
implicit none
double precision :: rpKC, rpAMP, phi, ev
integer :: i, j, k
ev = ((0.5*NX/PI)**2.0)*ENERV
rpAMP = sqrt(((0.11095279993106999d0*NV**2.5d0)*(NX*NY*NZ))/(ev**1.5d0))
rpKC = (1.666666666666666d0*ev)/NV
if (RANK .eq. 0) then
write (6, *) ' Imposing the highly non-equilibrium state...'
write (6, *) ' K-space amplitude = ', rpAMP
write (6, *) ' Maximum wavenumber = ', sqrt(rpKC)*DKSPACE
end if
!$OMP parallel do private (i,j,k) collapse(3)
do k = sz, ez
do j = sy, ey
do i = sx, ex
if (KX(i)**2 + KY(j)**2 + KZ(k)**2 <= rpKC*DKSPACE**2) then
call random_number(phi)
WS%FLUID(1)%GRID(i, j, k) = rpAMP*exp(2.0*PI*EYE*phi)
else
WS%FLUID(1)%GRID(i, j, k) = 0.0d0
end if
end do
end do
end do
!$OMP end parallel do
call fftw_mpi_execute_dft(fftw_backward_plan, WS%FLUID(1)%GRID, WS%FLUID(1)%GRID)
WS%FLUID(1)%GRID = WS%FLUID(1)%GRID/sqrt(dble(NX*NY*NZ))
end subroutine
subroutine insert_vortex_ring(field, x0, y0, z0, r0, dir, rot)
implicit none
integer :: i, j, k, dir, rot
type(fluid_field) :: field
double precision :: x0, y0, z0, r0
double precision, dimension(sx:ex, sy:ey, sz:ez) :: rr1, rr2, d1, d2
double precision, dimension(:, :), allocatable :: s
if (rot .eq. 0) then
allocate (s(sx:ex, sy:ey))
!$OMP parallel do private (i,j) collapse(2)
do j = sy, ey
do i = sx, ex
s(i, j) = sqrt((GX(i) - x0)**2 + (GY(j) - y0)**2)
end do
end do
!$OMP end parallel do
else if (rot .eq. 1) then
allocate (s(sx:ex, sz:ez))
!$OMP parallel do private (i,k) collapse(2)
do k = sz, ez
do i = sx, ex
s(i, k) = sqrt((GX(i) - x0)**2 + (GZ(k) - z0)**2)
end do
end do
!$OMP end parallel do
else if (rot .eq. 2) then
allocate (s(sy:ey, sz:ez))
!$OMP parallel do private (k,j) collapse(2)
do k = sz, ez
do j = sy, ey
s(j, k) = sqrt((GY(j) - y0)**2 + (GZ(k) - z0)**2)
end do
end do
!$OMP end parallel do
end if
!$OMP parallel do private (i,j,k) collapse(3)
do k = sz, ez
do j = sy, ey
do i = sx, ex
if (rot .eq. 0) then
d1(i, j, k) = sqrt((GZ(k) - z0)**2 + (s(i, j) + r0)**2)
d2(i, j, k) = sqrt((GZ(k) - z0)**2 + (s(i, j) - r0)**2)
else if (rot .eq. 1) then
d1(i, j, k) = sqrt((GY(j) - y0)**2 + (s(i, k) + r0)**2)
d2(i, j, k) = sqrt((GY(j) - y0)**2 + (s(i, k) - r0)**2)
else if (rot .eq. 2) then
d1(i, j, k) = sqrt((GX(i) - x0)**2 + (s(j, k) + r0)**2)
d2(i, j, k) = sqrt((GX(i) - x0)**2 + (s(j, k) - r0)**2)
end if
end do
end do
end do
!$OMP end parallel do
rr1 = sqrt(((0.3437d0 + 0.0572d0*d1**2))/(1.0d0 + (0.6666d0*d1**2) + (0.0572d0*d1**4)))
rr2 = sqrt(((0.3437d0 + 0.0572d0*d2**2))/(1.0d0 + (0.6666d0*d2**2) + (0.0572d0*d2**4)))
!$OMP parallel do private (i,j,k) collapse(3)
do k = sz, ez
do j = sy, ey
do i = sx, ex
if (rot .eq. 0) then
field%GRID(i, j, k) = field%GRID(i, j, k)*rr1(i, j, k)*(GZ(k) - z0 + dir*EYE*(s(i, j) + r0))* &
rr2(i, j, k)*(GZ(k) - z0 - dir*EYE*(s(i, j) - r0))
else if (rot .eq. 1) then
field%GRID(i, j, k) = field%GRID(i, j, k)*rr1(i, j, k)*(GY(j) - y0 + dir*EYE*(s(i, k) + r0))* &
rr2(i, j, k)*(GY(j) - y0 - dir*EYE*(s(i, k) - r0))
else if (rot .eq. 2) then
field%GRID(i, j, k) = field%GRID(i, j, k)*rr1(i, j, k)*(GX(i) - x0 + dir*EYE*(s(j, k) + r0))* &
rr2(i, j, k)*(GX(i) - x0 - dir*EYE*(s(j, k) - r0))
end if
end do
end do
end do
!$OMP end parallel do
deallocate (s)
end subroutine
end module