Skip to content

Commit 1980b32

Browse files
AmazingkivasAmazingkivas
committed
fix currents iterations
1 parent 4a85c45 commit 1980b32

2 files changed

Lines changed: 334 additions & 8 deletions

File tree

‎coarray/fdtd.F90‎

Lines changed: 31 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ program fdtd_coarray_optimized
2424
integer(int64) :: start_time, end_time, rate
2525
integer :: start_i, start_j, start_k, max_i, max_j, max_k
2626
real(8) :: elapsed_time
27-
real(8) :: Tx, Ty, Tz, TT, bnd
27+
real(8) :: Tx, Ty, Tz, TT, bnd, current_time
2828

2929
! Coarray initialization
3030
this_img = this_image()
@@ -36,6 +36,7 @@ program fdtd_coarray_optimized
3636
Tx = 4 * C
3737
Ty = 4 * C
3838
Tz = 4 * C
39+
current_time = min(int(TT / dt), num_iterations)
3940

4041
bnd = Ni / 2.0 * dx
4142

@@ -77,9 +78,7 @@ program fdtd_coarray_optimized
7778
end if
7879

7980
! Main FDTD loop
80-
do t = 1, num_iterations
81-
if (this_img == 1) print '(A, I3)', 'Iteration ', t
82-
81+
do t = 1, current_time
8382
if (this_img == 1) then
8483
call init_currents(t)
8584
end if
@@ -103,6 +102,31 @@ program fdtd_coarray_optimized
103102
sync all
104103

105104
end do
105+
106+
Jx = 0.0
107+
Jy = 0.0
108+
Jz = 0.0
109+
sync all
110+
111+
do t = current_time + 1, num_iterations
112+
call update_B_field()
113+
sync all
114+
115+
pred_Bx(:,:) = Bx(:, :, k_local)[pred_img]
116+
pred_By(:,:) = By(:, :, k_local)[pred_img]
117+
sync all
118+
119+
call update_E_field()
120+
sync all
121+
122+
next_Ex(:,:) = Ex(:, :, 1)[next_img]
123+
next_Ey(:,:) = Ey(:, :, 1)[next_img]
124+
sync all
125+
126+
call update_B_field()
127+
sync all
128+
129+
end do
106130

107131
sync all
108132
! Final timing and cleanup
@@ -267,11 +291,10 @@ end subroutine update_E_field
267291
!===============================================================
268292
subroutine print_full_E_slice()
269293
integer :: img, ik, ij
270-
do ij = Nj/2-5, Nj/2+5
271-
do ik = Ni/2-5, Ni/2+5
294+
do ij = Nj/2-4, Nj/2+5
295+
do ik = Ni/2-4, Ni/2+5
272296
img = ceiling(real(Nk/2+1) / real(k_local))
273-
write(*, '(F12.6)', advance='no') Ex(ik,ij,Nk/2+1 - (img - 1) * k_local)[img]
274-
! write(*, '(F12.6)', advance='no') next_Ex(ik, ij)
297+
write(*, '(F12.5)', advance='no') Ex(ik,ij,Nk/2+1 - (img - 1) * k_local)[img]
275298
end do
276299
print *
277300
end do

‎coarray/fdtd.f90‎

Lines changed: 303 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,303 @@
1+
program fdtd_coarray_optimized
2+
use iso_fortran_env
3+
implicit none
4+
5+
! Parameters
6+
integer, parameter :: Ni = 32, Nj = 32, Nk = 32
7+
integer, parameter :: num_iterations = 100
8+
real(8), parameter :: C = 3e10, PI = 3.14159265358
9+
real(8), parameter :: dx = C, dy = C, dz = C, dt = 0.2
10+
real(8), parameter :: coef_B_dx = C * dt / (2 * dx), coef_B_dy = C * dt / (2 * dy), coef_B_dz = C * dt / (2 * dz)
11+
real(8), parameter :: coef_E_dx = C * dt / dx, coef_E_dy = C * dt / dy, coef_E_dz = C * dt / dz
12+
real(8), parameter :: coef_J = 4*PI*dt
13+
14+
! Field arrays (optimized layout for vectorization)
15+
real(8), allocatable :: Ex(:,:,:)[:], Ey(:,:,:)[:], Ez(:,:,:)[:]
16+
real(8), allocatable :: Bx(:,:,:)[:], By(:,:,:)[:], Bz(:,:,:)[:]
17+
real(8), allocatable :: Jx(:,:,:)[:], Jy(:,:,:)[:], Jz(:,:,:)[:]
18+
real(8), allocatable :: next_Ex(:,:)[:], next_Ey(:,:)[:], pred_Bx(:,:)[:], pred_By(:,:)[:]
19+
20+
! Local variables
21+
integer :: this_img, total_imgs, next_img, pred_img
22+
integer :: k_start, k_end, k_local
23+
integer :: idx, t
24+
integer(int64) :: start_time, end_time, rate
25+
integer :: start_i, start_j, start_k, max_i, max_j, max_k
26+
real(8) :: elapsed_time
27+
real(8) :: Tx, Ty, Tz, TT, bnd, current_time
28+
29+
! Coarray initialization
30+
this_img = this_image()
31+
total_imgs = num_images()
32+
next_img = merge(this_image() + 1, 1, this_image() < num_images())
33+
pred_img = merge(this_image() - 1, num_images(), this_image() > 1)
34+
35+
TT = 8.0
36+
Tx = 4 * C
37+
Ty = 4 * C
38+
Tz = 4 * C
39+
current_time = min(int(TT / dt), num_iterations)
40+
41+
bnd = Ni / 2.0 * dx
42+
43+
start_i = floor((-Tx/4.0 + bnd) / dx) + 1
44+
start_j = floor((-Ty/4.0 + bnd) / dy) + 1
45+
start_k = floor((-Tz/4.0 + bnd) / dz) + 1
46+
47+
max_i = floor((Tx/4.0 + bnd) / dx) + 1
48+
max_j = floor((Ty/4.0 + bnd) / dy) + 1
49+
max_k = floor((Tz/4.0 + bnd) / dz) + 1
50+
51+
! Domain decomposition
52+
call init_decomposition()
53+
54+
! Allocate arrays with optimized layout
55+
allocate(Ex(Ni, Nj, k_local)[*])
56+
allocate(Ey(Ni, Nj, k_local)[*])
57+
allocate(Ez(Ni, Nj, k_local)[*])
58+
allocate(Bx(Ni, Nj, k_local)[*])
59+
allocate(By(Ni, Nj, k_local)[*])
60+
allocate(Bz(Ni, Nj, k_local)[*])
61+
allocate(Jx(Ni, Nj, k_local)[*])
62+
allocate(Jy(Ni, Nj, k_local)[*])
63+
allocate(Jz(Ni, Nj, k_local)[*])
64+
65+
allocate(next_Ex(Ni, Nj)[*])
66+
allocate(next_Ey(Ni, Nj)[*])
67+
allocate(pred_Bx(Ni, Nj)[*])
68+
allocate(pred_By(Ni, Nj)[*])
69+
70+
! Initialize fields
71+
call init_fields()
72+
73+
! Timing and output
74+
if (this_img == 1) then
75+
print '(A, I0)', 'Running on ', total_imgs, ' images'
76+
call system_clock(count_rate=rate)
77+
call system_clock(count=start_time)
78+
end if
79+
80+
! Main FDTD loop
81+
do t = 1, current_time
82+
if (this_img == 1) then
83+
call init_currents(t)
84+
end if
85+
sync all
86+
87+
call update_B_field()
88+
sync all
89+
90+
pred_Bx(:,:) = Bx(:, :, k_local)[pred_img]
91+
pred_By(:,:) = By(:, :, k_local)[pred_img]
92+
sync all
93+
94+
call update_E_field()
95+
sync all
96+
97+
next_Ex(:,:) = Ex(:, :, 1)[next_img]
98+
next_Ey(:,:) = Ey(:, :, 1)[next_img]
99+
sync all
100+
101+
call update_B_field()
102+
sync all
103+
104+
end do
105+
106+
Jx = 0.0
107+
Jy = 0.0
108+
Jz = 0.0
109+
sync all
110+
111+
do t = current_time + 1, num_iterations
112+
call update_B_field()
113+
sync all
114+
115+
pred_Bx(:,:) = Bx(:, :, k_local)[pred_img]
116+
pred_By(:,:) = By(:, :, k_local)[pred_img]
117+
sync all
118+
119+
call update_E_field()
120+
sync all
121+
122+
next_Ex(:,:) = Ex(:, :, 1)[next_img]
123+
next_Ey(:,:) = Ey(:, :, 1)[next_img]
124+
sync all
125+
126+
call update_B_field()
127+
sync all
128+
129+
end do
130+
131+
sync all
132+
! Final timing and cleanup
133+
if (this_img == 1) then
134+
call system_clock(count=end_time)
135+
elapsed_time = real(end_time - start_time)/real(rate)
136+
print '(A, F0.2, A)', 'Total execution time: ', elapsed_time, ' seconds'
137+
end if
138+
139+
sync all
140+
if ((this_image() == 1)) then !.and. (Nk <= 16)) then
141+
call print_full_E_slice()
142+
end if
143+
144+
deallocate(Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz)
145+
deallocate(next_Ex, next_Ey, pred_Bx, pred_By)
146+
147+
contains
148+
!===============================================================
149+
subroutine init_decomposition
150+
integer :: remainder, offset
151+
152+
remainder = mod(Nk, total_imgs)
153+
k_local = Nk / total_imgs
154+
155+
if (this_img <= remainder) then
156+
k_local = k_local + 1
157+
offset = (this_img - 1)*k_local
158+
else
159+
offset = remainder*(k_local + 1) + (this_img - remainder - 1)*k_local
160+
end if
161+
162+
k_start = offset + 1
163+
k_end = offset + k_local
164+
end subroutine init_decomposition
165+
166+
!===============================================================
167+
subroutine init_fields
168+
Ex = 0.0
169+
Ey = 0.0
170+
Ez = 0.0
171+
next_Ex = 0.0
172+
next_Ey = 0.0
173+
pred_Bx = 0.0
174+
pred_By = 0.0
175+
Bx = 0.0
176+
By = 0.0
177+
Bz = 0.0
178+
Jx = 0.0
179+
Jy = 0.0
180+
Jz = 0.0
181+
end subroutine init_fields
182+
183+
!===============================================================
184+
subroutine init_currents(this_t)
185+
integer, intent(in) :: this_t
186+
integer :: i, j, k
187+
integer :: this_cur_img
188+
real(8) :: value
189+
190+
do k = start_k, max_k
191+
do j = start_j, max_j
192+
do i = start_i, max_i
193+
this_cur_img = ceiling(real(k) / real(k_local))
194+
value = (sin(2.0 * PI * this_t * dt / TT)) &
195+
* (cos(2.0 * PI * (i-1) * dx / Tx)**2) &
196+
* (cos(2.0 * PI * (j-1) * dy / Ty)**2) &
197+
* (cos(2.0 * PI * (k-1) * dz / Tz)**2)
198+
Jx(i,j,k - (this_cur_img - 1) * k_local)[this_cur_img] = value
199+
Jy(i,j,k - (this_cur_img - 1) * k_local)[this_cur_img] = value
200+
Jz(i,j,k - (this_cur_img - 1) * k_local)[this_cur_img] = value
201+
end do
202+
end do
203+
end do
204+
end subroutine init_currents
205+
206+
!===============================================================
207+
subroutine update_B_field
208+
integer :: ip, jp, kp
209+
integer :: i, j, k
210+
211+
do k = 1, k_local - 1
212+
do j = 1, Nj
213+
jp = merge(j+1, 1, j < Nj)
214+
do i = 1, Ni
215+
ip = merge(i+1, 1, i < Ni)
216+
217+
Bx(i,j,k) = Bx(i,j,k) + coef_B_dz * (Ey(i,j,k+1) - Ey(i,j,k)) - &
218+
coef_B_dy * (Ez(i,jp,k) - Ez(i,j,k))
219+
220+
By(i,j,k) = By(i,j,k) + coef_B_dx * (Ez(ip,j,k) - Ez(i,j,k)) - &
221+
coef_B_dz * (Ex(i,j,k+1) - Ex(i,j,k))
222+
223+
Bz(i,j,k) = Bz(i,j,k) + coef_B_dy * (Ex(i,jp,k) - Ex(i,j,k)) - &
224+
coef_B_dx * (Ey(ip,j,k) - Ey(i,j,k))
225+
end do
226+
end do
227+
end do
228+
do j = 1, Nj
229+
jp = merge(j+1, 1, j < Nj)
230+
do i = 1, Ni
231+
ip = merge(i+1, 1, i < Ni)
232+
233+
Bx(i,j,k_local) = Bx(i,j,k_local) + coef_B_dz * (next_Ey(i,j) - Ey(i,j,k_local)) - &
234+
coef_B_dy * (Ez(i,jp,k_local) - Ez(i,j,k_local))
235+
236+
By(i,j,k_local) = By(i,j,k_local) + coef_B_dx * (Ez(ip,j,k_local) - Ez(i,j,k)) - &
237+
coef_B_dz * (next_Ex(i,j) - Ex(i,j,k_local))
238+
239+
Bz(i,j,k_local) = Bz(i,j,k_local) + coef_B_dy * (Ex(i,jp,k_local) - Ex(i,j,k_local)) - &
240+
coef_B_dx * (Ey(ip,j,k_local) - Ey(i,j,k_local))
241+
end do
242+
end do
243+
244+
end subroutine update_B_field
245+
246+
!===============================================================
247+
subroutine update_E_field
248+
integer :: im, jm, km
249+
integer :: i, j, k
250+
251+
do j = 1, Nj
252+
jm = merge(j-1, Nj, j > 1)
253+
do i = 1, Ni
254+
im = merge(i-1, Ni, i > 1)
255+
256+
Ex(i,j,1) = Ex(i,j,1) - coef_J * Jx(i,j,1) + &
257+
coef_E_dy * (Bz(i,j,1) - Bz(i,jm,1)) - &
258+
coef_E_dz * (By(i,j,1) - pred_By(i,j))
259+
260+
Ey(i,j,1) = Ey(i,j,1) - coef_J * Jy(i,j,1) + &
261+
coef_E_dz * (Bx(i,j,1) - pred_Bx(i,j)) - &
262+
coef_E_dx * (Bz(i,j,1) - Bz(im,j,1))
263+
264+
Ez(i,j,1) = Ez(i,j,1) - coef_J * Jz(i,j,1) + &
265+
coef_E_dx * (By(i,j,1) - By(im,j,1)) - &
266+
coef_E_dy * (Bx(i,j,1) - Bx(i,jm,1))
267+
end do
268+
end do
269+
do k = 2, k_local
270+
do j = 1, Nj
271+
jm = merge(j-1, Nj, j > 1)
272+
do i = 1, Ni
273+
im = merge(i-1, Ni, i > 1)
274+
275+
Ex(i,j,k) = Ex(i,j,k) - coef_J * Jx(i,j,k) + &
276+
coef_E_dy * (Bz(i,j,k) - Bz(i,jm,k)) - &
277+
coef_E_dz * (By(i,j,k) - By(i,j,k-1))
278+
279+
Ey(i,j,k) = Ey(i,j,k) - coef_J * Jy(i,j,k) + &
280+
coef_E_dz * (Bx(i,j,k) - Bx(i,j,k-1)) - &
281+
coef_E_dx * (Bz(i,j,k) - Bz(im,j,k))
282+
283+
Ez(i,j,k) = Ez(i,j,k) - coef_J * Jz(i,j,k) + &
284+
coef_E_dx * (By(i,j,k) - By(im,j,k)) - &
285+
coef_E_dy * (Bx(i,j,k) - Bx(i,jm,k))
286+
end do
287+
end do
288+
end do
289+
end subroutine update_E_field
290+
291+
!===============================================================
292+
subroutine print_full_E_slice()
293+
integer :: img, ik, ij
294+
do ij = Nj/2-4, Nj/2+5
295+
do ik = Ni/2-4, Ni/2+5
296+
img = ceiling(real(Nk/2+1) / real(k_local))
297+
write(*, '(F12.5)', advance='no') Ex(ik,ij,Nk/2+1 - (img - 1) * k_local)[img]
298+
end do
299+
print *
300+
end do
301+
302+
end subroutine print_full_E_slice
303+
end program fdtd_coarray_optimized

0 commit comments

Comments
 (0)