Skip to content

Commit 8687aaf

Browse files
HarryMcleanHarryMclean
committed
Enable threading
1 parent b2b58f5 commit 8687aaf

3 files changed

Lines changed: 23 additions & 73 deletions

File tree

‎src/heatflow.f90‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,12 @@ program HEATFLOW_V0_3
3434
real(real12) :: cpustart, cpuend, cpustart2, progress
3535
integer(int12) :: itime
3636

37+
!-------------------------------------------------------------!
38+
! Initialize PETSc FIRST (before any other operations) !
39+
!-------------------------------------------------------------!
40+
CALL petsc_init()
41+
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
42+
3743
!-------------------------------------------------------------!
3844
! calculate the time to run full simulation !
3945
!-------------------------------------------------------------!
@@ -65,7 +71,6 @@ program HEATFLOW_V0_3
6571
!-------------------------------------------------------------!
6672
! run simulation for 'ntime' time steps !
6773
!-------------------------------------------------------------!
68-
CALL petsc_init()
6974

7075
do itime=1,ntime
7176

@@ -99,7 +104,7 @@ program HEATFLOW_V0_3
99104
!-------------------------------------------------------------!
100105
! calculate end time and print to user !
101106
!-------------------------------------------------------------!
102-
CALL cpu_time(cpuend)
107+
CALL cpu_time(cpuend)
103108
write(*,'(A,F12.6)') ' time=', cpuend-cpustart
104109
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
105110

‎src/heatflow/mod_petsc_solver.f90‎

Lines changed: 14 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -60,42 +60,42 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
6060
real(8) :: rnorm
6161
logical :: rebuild_needed
6262

63-
write(*,'(A,I0)') ' [DEBUG] Entered solve_petsc_csr, n=', n
64-
call flush(6)
63+
! write(*,'(A,I0)') ' [DEBUG] Entered solve_petsc_csr, n=', n
64+
! call flush(6)
6565

6666
if (size(ia) /= n+1) stop 'solve_petsc_csr: ia size mismatch'
6767
if (size(b) /= n .or. size(x) /= n) stop 'solve_petsc_csr: vector size mismatch'
6868

69-
write(*,'(A)') ' [DEBUG] Size checks passed'
70-
call flush(6)
69+
! write(*,'(A)') ' [DEBUG] Size checks passed'
70+
! call flush(6)
7171

7272
! Determine if we need to rebuild the matrix structure
7373
rebuild_needed = .false.
7474
if (.not. initialized) rebuild_needed = .true.
7575
if (n /= n_saved) rebuild_needed = .true.
7676

77-
write(*,'(A,L1)') ' [DEBUG] rebuild_needed=', rebuild_needed
78-
call flush(6)
77+
! write(*,'(A,L1)') ' [DEBUG] rebuild_needed=', rebuild_needed
78+
! call flush(6)
7979

8080
! Create PETSc objects on first call or if size changed
8181
if (rebuild_needed) then
82-
write(*,'(A)') ' [DEBUG] Starting PETSc object creation...'
83-
call flush(6)
82+
! write(*,'(A)') ' [DEBUG] Starting PETSc object creation...'
83+
! call flush(6)
8484

8585
! Clean up old objects if they exist
8686
if (initialized) call petsc_cleanup()
8787

88-
write(*,'(A)') ' [DEBUG] Preallocating matrix...'
89-
call flush(6)
88+
! write(*,'(A)') ' [DEBUG] Preallocating matrix...'
89+
! call flush(6)
9090

9191
! Preallocate matrix with exact nonzeros per row (saves memory)
9292
allocate(d_nnz(n))
9393
do i = 1, n
9494
d_nnz(i) = ia(i+1) - ia(i)
9595
end do
9696

97-
write(*,'(A,I0,A,I0)') ' [DEBUG] Creating matrix: n=', n, ', max_nnz/row=', maxval(d_nnz)
98-
call flush(6)
97+
! write(*,'(A,I0,A,I0)') ' [DEBUG] Creating matrix: n=', n, ', max_nnz/row=', maxval(d_nnz)
98+
! call flush(6)
9999

100100
! Create matrix with exact preallocation (most memory-efficient)
101101
call MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 0, d_nnz, A_saved, ierr)
@@ -106,18 +106,12 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
106106
stop
107107
end if
108108

109-
write(*,'(A)') ' [DEBUG] Matrix created successfully'
110-
call flush(6)
111-
112109
deallocate(d_nnz)
113110

114111
! Create persistent vectors
115112
call VecCreateSeq(PETSC_COMM_SELF, n, bb_saved, ierr)
116113
call VecCreateSeq(PETSC_COMM_SELF, n, xx_saved, ierr)
117114

118-
write(*,'(A)') ' [DEBUG] Vectors created'
119-
call flush(6)
120-
121115
! Create and configure KSP solver (persistent across timesteps)
122116
call KSPCreate(PETSC_COMM_SELF, ksp_saved, ierr)
123117
call KSPSetOperators(ksp_saved, A_saved, A_saved, ierr)
@@ -134,8 +128,6 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
134128
end if
135129

136130
! Update matrix values (always needed each timestep)
137-
write(*,'(A)') ' Updating PETSc matrix values...'
138-
call flush(6)
139131
call MatZeroEntries(A_saved, ierr)
140132
do i = 1, n
141133
row_nz = ia(i+1) - ia(i)
@@ -151,19 +143,13 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
151143
deallocate(cols0, vals)
152144
end if
153145
end do
154-
write(*,'(A)') ' Matrix assembly beginning...'
155-
call flush(6)
156146
call MatAssemblyBegin(A_saved, MAT_FINAL_ASSEMBLY, ierr)
157147
call MatAssemblyEnd(A_saved, MAT_FINAL_ASSEMBLY, ierr)
158-
write(*,'(A)') ' Matrix assembly complete.'
159-
call flush(6)
160148

161149
! Optional: Verify matrix assembly (uncomment for debugging)
162150
! call MatView(A_saved, PETSC_VIEWER_STDOUT_SELF, ierr)
163151

164152
! Update RHS vector in batches to avoid memory issues with huge systems
165-
write(*,'(A)') ' Updating RHS vector...'
166-
call flush(6)
167153
block
168154
integer, parameter :: VEC_CHUNK = 1000000
169155
integer :: vec_start, vec_end, vec_len, k
@@ -177,21 +163,12 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
177163
idx_vec = [(vec_start + k - 2, k=1,vec_len)] ! 0-based indices
178164

179165
call VecSetValues(bb_saved, vec_len, idx_vec, b(vec_start:vec_end), INSERT_VALUES, ierr)
180-
if (ierr /= 0) then
181-
write(0,*) "ERROR: VecSetValues(bb) failed at", vec_start, "ierr=", ierr
182-
stop
183-
end if
184-
185166
deallocate(idx_vec)
186167
end do
187168
end block
188169
call VecAssemblyBegin(bb_saved,ierr); call VecAssemblyEnd(bb_saved,ierr)
189-
write(*,'(A)') ' RHS vector complete.'
190-
call flush(6)
191170

192171
! Update initial guess in batches
193-
write(*,'(A)') ' Updating initial guess vector...'
194-
call flush(6)
195172
block
196173
integer, parameter :: VEC_CHUNK = 1000000
197174
integer :: vec_start, vec_end, vec_len, k
@@ -205,49 +182,23 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
205182
idx_vec = [(vec_start + k - 2, k=1,vec_len)] ! 0-based indices
206183

207184
call VecSetValues(xx_saved, vec_len, idx_vec, x(vec_start:vec_end), INSERT_VALUES, ierr)
208-
if (ierr /= 0) then
209-
write(0,*) "ERROR: VecSetValues(xx) failed at", vec_start, "ierr=", ierr
210-
stop
211-
end if
212-
213185
deallocate(idx_vec)
214186
end do
215187
end block
216188
call VecAssemblyBegin(xx_saved,ierr); call VecAssemblyEnd(xx_saved,ierr)
217189

218-
! Solve the system (reusing persistent KSP)
219-
write(*,'(A,I0,A)') ' Solving linear system with n=', n, ' unknowns...'
220-
call flush(6)
190+
! Solve the system
221191
call KSPSolve(ksp_saved, bb_saved, xx_saved, ierr)
222192

223-
write(*,'(A)') ' KSPSolve completed, checking status...'
224-
call flush(6)
225-
226193
if (ierr /= 0) then
227194
write(0,*) "ERROR: KSPSolve failed with error code:", ierr
228195
stop
229196
end if
230197

231-
! Check convergence - Note: KSPConvergedReason type changed in PETSc 3.24
232-
! Simplified error checking without explicit reason query
233198
call KSPGetIterationNumber(ksp_saved, its, ierr)
234199
call KSPGetResidualNorm(ksp_saved, rnorm, ierr)
235-
236-
! Report convergence status (commented out by default for performance)
237-
! Uncomment the next line to see convergence info every solve:
238-
! write(*,'(A,I0,A,ES12.5)') ' PETSc: iterations=', its, ', residual=', rnorm
239-
240-
! Basic divergence check via error code from solve
241-
if (ierr /= 0) then
242-
write(0,*) "WARNING: PETSc solver returned non-zero error code!"
243-
write(0,*) " Error code:", ierr
244-
write(0,*) " Iterations:", its
245-
write(0,*) " Residual norm:", rnorm
246-
! Don't stop - let the main code detect NaNs if needed
247-
end if
248200

249201
! Extract solution vector using batched VecGetValues
250-
! Process in chunks for better performance on large systems
251202
block
252203
integer, parameter :: CHUNK_SIZE = 100000
253204
PetscInt, allocatable :: idx_batch(:)
@@ -260,19 +211,11 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
260211

261212
allocate(idx_batch(chunk_len), val_batch(chunk_len))
262213

263-
! Build index array (0-based for PETSc)
264214
do j = 1, chunk_len
265-
idx_batch(j) = i_start + j - 2 ! -1 for 0-based, -1 more for offset
215+
idx_batch(j) = i_start + j - 2
266216
end do
267217

268-
! Get chunk of values
269218
call VecGetValues(xx_saved, chunk_len, idx_batch, val_batch, ierr)
270-
if (ierr /= 0) then
271-
write(0,*) "ERROR: VecGetValues failed at chunk starting", i_start, "ierr=", ierr
272-
stop
273-
end if
274-
275-
! Copy to output array
276219
x(i_start:i_end) = val_batch(1:chunk_len)
277220

278221
deallocate(idx_batch, val_batch)

‎src/heatflow/mod_setup.f90‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,8 @@ module setup
2626
use materials, only: material
2727
implicit none
2828

29+
public :: set_global_variables
30+
2931
contains
3032

3133

0 commit comments

Comments
 (0)