@@ -3,8 +3,8 @@ program fdtd_coarray_optimized
33 implicit none
44
55 ! Parameters
6- integer , parameter :: Ni = 512 , Nj = 512 , Nk = 512
7- integer , parameter :: num_iterations = 25
6+ integer , parameter :: Ni = 32 , Nj = 32 , Nk = 32
7+ integer , parameter :: num_iterations = 100
88 real (8 ), parameter :: C = 3e10 , PI = 3.14159265358
99 real (8 ), parameter :: dx = C, dy = C, dz = C, dt = 0.2
1010 real (8 ), parameter :: coef_B_dx = C * dt / (2 * dx), coef_B_dy = C * dt / (2 * dy), coef_B_dz = C * dt / (2 * dz)
@@ -14,7 +14,7 @@ program fdtd_coarray_optimized
1414 ! Field arrays (optimized layout for vectorization)
1515 real (8 ), allocatable :: Ex(:,:,:)[:], Ey(:,:,:)[:], Ez(:,:,:)[:]
1616 real (8 ), allocatable :: Bx(:,:,:)[:], By(:,:,:)[:], Bz(:,:,:)[:]
17- real (8 ), allocatable :: Jx(:,:,:)[:]
17+ real (8 ), allocatable :: Jx(:,:,:)[:], Jy(:,:,:)[:], Jz(:,:,:)[:]
1818 real (8 ), allocatable :: next_Ex(:,:)[:], next_Ey(:,:)[:], pred_Bx(:,:)[:], pred_By(:,:)[:]
1919
2020 ! Local variables
@@ -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
@@ -58,7 +59,9 @@ program fdtd_coarray_optimized
5859 allocate (By(Ni, Nj, k_local)[* ])
5960 allocate (Bz(Ni, Nj, k_local)[* ])
6061 allocate (Jx(Ni, Nj, k_local)[* ])
61-
62+ allocate (Jy(Ni, Nj, k_local)[* ])
63+ allocate (Jz(Ni, Nj, k_local)[* ])
64+
6265 allocate (next_Ex(Ni, Nj)[* ])
6366 allocate (next_Ey(Ni, Nj)[* ])
6467 allocate (pred_Bx(Ni, Nj)[* ])
@@ -75,9 +78,7 @@ program fdtd_coarray_optimized
7578 end if
7679
7780 ! Main FDTD loop
78- do t = 1 , num_iterations
79- if (this_img == 1 ) print ' (A, I3)' , ' Iteration ' , t
80-
81+ do t = 1 , current_time
8182 if (this_img == 1 ) then
8283 call init_currents(t)
8384 end if
@@ -101,6 +102,31 @@ program fdtd_coarray_optimized
101102 sync all
102103
103104 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
104130
105131 sync all
106132 ! Final timing and cleanup
@@ -115,7 +141,7 @@ program fdtd_coarray_optimized
115141 call print_full_E_slice()
116142 end if
117143
118- deallocate (Ex, Ey, Ez, Bx, By, Bz, Jx)
144+ deallocate (Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz )
119145 deallocate (next_Ex, next_Ey, pred_Bx, pred_By)
120146
121147contains
@@ -150,23 +176,28 @@ subroutine init_fields
150176 By = 0.0
151177 Bz = 0.0
152178 Jx = 0.0
179+ Jy = 0.0
180+ Jz = 0.0
153181 end subroutine init_fields
154182
155183 ! ===============================================================
156184 subroutine init_currents (this_t )
157185 integer , intent (in ) :: this_t
158186 integer :: i, j, k
159187 integer :: this_cur_img
188+ real (8 ) :: value
160189
161190 do k = start_k, max_k
162191 do j = start_j, max_j
163192 do i = start_i, max_i
164193 this_cur_img = ceiling (real (k) / real (k_local))
165- Jx(i,j,k - (this_cur_img - 1 ) * k_local)[this_cur_img] = &
166- (sin (2.0 * PI * this_t * dt / TT)) &
167- * (cos (2.0 * PI * (i-1 ) * dx / Tx)** 2 ) &
168- * (cos (2.0 * PI * (j-1 ) * dy / Ty)** 2 ) &
169- * (cos (2.0 * PI * (k-1 ) * dz / Tz)** 2 )
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
170201 end do
171202 end do
172203 end do
@@ -226,11 +257,11 @@ subroutine update_E_field
226257 coef_E_dy * (Bz(i,j,1 ) - Bz(i,jm,1 )) - &
227258 coef_E_dz * (By(i,j,1 ) - pred_By(i,j))
228259
229- Ey(i,j,1 ) = Ey(i,j,1 ) - coef_J * Jx (i,j,1 ) + &
260+ Ey(i,j,1 ) = Ey(i,j,1 ) - coef_J * Jy (i,j,1 ) + &
230261 coef_E_dz * (Bx(i,j,1 ) - pred_Bx(i,j)) - &
231262 coef_E_dx * (Bz(i,j,1 ) - Bz(im,j,1 ))
232263
233- Ez(i,j,1 ) = Ez(i,j,1 ) - coef_J * Jx (i,j,1 ) + &
264+ Ez(i,j,1 ) = Ez(i,j,1 ) - coef_J * Jz (i,j,1 ) + &
234265 coef_E_dx * (By(i,j,1 ) - By(im,j,1 )) - &
235266 coef_E_dy * (Bx(i,j,1 ) - Bx(i,jm,1 ))
236267 end do
@@ -245,11 +276,11 @@ subroutine update_E_field
245276 coef_E_dy * (Bz(i,j,k) - Bz(i,jm,k)) - &
246277 coef_E_dz * (By(i,j,k) - By(i,j,k-1 ))
247278
248- Ey(i,j,k) = Ey(i,j,k) - coef_J * Jx (i,j,k) + &
279+ Ey(i,j,k) = Ey(i,j,k) - coef_J * Jy (i,j,k) + &
249280 coef_E_dz * (Bx(i,j,k) - Bx(i,j,k-1 )) - &
250281 coef_E_dx * (Bz(i,j,k) - Bz(im,j,k))
251282
252- Ez(i,j,k) = Ez(i,j,k) - coef_J * Jx (i,j,k) + &
283+ Ez(i,j,k) = Ez(i,j,k) - coef_J * Jz (i,j,k) + &
253284 coef_E_dx * (By(i,j,k) - By(im,j,k)) - &
254285 coef_E_dy * (Bx(i,j,k) - Bx(i,jm,k))
255286 end do
@@ -260,10 +291,10 @@ end subroutine update_E_field
260291 ! ===============================================================
261292subroutine print_full_E_slice ()
262293 integer :: img, ik, ij
263- do ik = Ni / 2-5 , Ni / 2+5
264- do ij = Nj / 2-5 , Nj / 2+5
294+ do ij = Nj / 2-4 , Nj / 2+5
295+ do ik = Ni / 2-4 , Ni / 2+5
265296 img = ceiling (real (Nk/ 2+1 ) / real (k_local))
266- write (* , ' (F12.6 )' , advance= ' no' ) Ex(ik,ij,Nk/ 2+1 - (img - 1 ) * k_local)[img]
297+ write (* , ' (F12.5 )' , advance= ' no' ) Ex(ik,ij,Nk/ 2+1 - (img - 1 ) * k_local)[img]
267298 end do
268299 print *
269300 end do
0 commit comments