Skip to content

Commit b8d7725

Browse files
fd9034fd9034
committed
Fix HDF5 output dimension ordering for Python compatibility (nt,nz,ny,nx). Add debug instrumentation.
1 parent bed4f86 commit b8d7725

11 files changed

Lines changed: 326 additions & 71 deletions

‎.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,3 +23,4 @@ build.sh
2323
test/test_run2/restart/TempDis.dat
2424
.gitignore
2525
test/test_run2/restart/TempDisTPD.dat
26+
tmp/

‎Makefile‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -57,8 +57,15 @@ endif
5757
# BLAS/LAPACK: Use Apple Accelerate on macOS, OpenBLAS on Linux
5858
UNAME_S := $(shell uname -s)
5959
ifeq ($(UNAME_S),Darwin)
60+
# macOS SDK sysroot (fixes 'library System not found' with Homebrew gfortran)
61+
MACOS_SDK := $(shell xcrun --show-sdk-path 2>/dev/null)
62+
ifneq ($(MACOS_SDK),)
63+
SYSROOT_FLAGS := -L$(MACOS_SDK)/usr/lib -F$(MACOS_SDK)/System/Library/Frameworks
64+
else
65+
SYSROOT_FLAGS :=
66+
endif
6067
# Apple Accelerate framework - optimized for Apple Silicon
61-
BLAS_FLAGS := -framework Accelerate -lgomp -lpthread -lm
68+
BLAS_FLAGS := $(SYSROOT_FLAGS) -framework Accelerate -lgomp -lpthread -lm
6269
BLAS_NOTE := (Apple Accelerate)
6370
else
6471
# Linux: Use OpenBLAS for multi-threaded BLAS/LAPACK

‎src/.DS_Store‎

6 KB
Binary file not shown.

‎src/heatflow/mod_boundary.f90‎

Lines changed: 34 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -79,48 +79,48 @@ subroutine boundary(B)
7979
kappa = grid(ix, iy, iz)%kappa
8080

8181
if (.not. Periodicx) then
82-
if (CylindricalGrid) then
83-
!-------------------------------------------------------
84-
! Cylindrical radial boundaries:
85-
! ix=1: center symmetry -> zero flux, no bath term
86-
! ix=nx: outer radius -> use kappaBoundNr, T_BathNr
87-
!-------------------------------------------------------
88-
if (ix .eq. 1) then
89-
! Symmetry at r=0: no boundary flux contribution
90-
! (B(I) unchanged, zero flux)
91-
end if
92-
if (ix .eq. nx) then
93-
kappaHarm = (2*kappa*kappaBoundNr/(kappa+kappaBoundNr)) / &
94-
(grid(ix, iy, iz)%Length(1)**2)
95-
if (kappa .ne. kappaBoundNr) kappaHarm = kappaHarm*BR
96-
! Apply cylindrical area correction: r_outer / r_center
97-
r_center = real(ix,real12) - 0.5_real12
98-
r_iface = real(ix,real12)
99-
kappaHarm = kappaHarm * r_iface / r_center
100-
B(I) = B(I) + (kappaHarm) * T_BathNr
101-
end if
102-
else
103-
if (ix .eq. 1) then
104-
if (CG_x_m) then
82+
clyBC:if (CylindricalGrid) then
83+
!-------------------------------------------------------
84+
! Cylindrical radial boundaries:
85+
! ix=1: center symmetry -> zero flux, no bath term
86+
! ix=nx: outer radius -> use kappaBoundNr, T_BathNr
87+
!-------------------------------------------------------
88+
if (ix .eq. 1) then
89+
! Symmetry at r=0: no boundary flux contribution
90+
B(I) = 0.0_real12
91+
end if
92+
if (ix .eq. nx) then
93+
kappaHarm = (2*kappa*kappaBoundNr/(kappa+kappaBoundNr)) / &
94+
(grid(ix, iy, iz)%Length(1)**2)
95+
if (kappa .ne. kappaBoundNr) kappaHarm = kappaHarm*BR
96+
! Apply cylindrical area correction: r_outer / r_center
97+
r_center = real(ix,real12) - 0.5_real12
98+
r_iface = real(ix,real12)
99+
kappaHarm = kappaHarm * r_iface / r_center
100+
B(I) = B(I) + (kappaHarm) * T_BathNr
101+
end if
102+
else
103+
if (ix .eq. 1) then
104+
if (CG_x_m) then
105105
B(I) = x1_power_dens
106-
else
106+
else
107107
kappaHarm = (2*kappa*kappaBoundx1/(kappa+kappaBoundx1)) / &
108-
(grid(ix, iy, iz)%Length(1)**2)
108+
(grid(ix, iy, iz)%Length(1)**2)
109109
if (kappa .ne. kappaBoundx1) kappaHarm = kappaHarm*BR
110110
B(I) = B(I) + (kappaHarm) * T_Bathx1 !+ boundray_term_vel(1_int12,iy,iz,T_Bathx1)
111-
end if
112-
end if
113-
if (ix .eq. nx) then
114-
if (CG_x_p) then
111+
end if
112+
end if
113+
if (ix .eq. nx) then
114+
if (CG_x_p) then
115115
B(I) = xn_power_dens
116-
else
116+
else
117117
kappaHarm = (2*kappa*kappaBoundNx/(kappa+kappaBoundNx)) / &
118-
(grid(ix, iy, iz)%Length(1)**2)
118+
(grid(ix, iy, iz)%Length(1)**2)
119119
if (kappa .ne. kappaBoundNx) kappaHarm = kappaHarm*BR
120120
B(I) = B(I) + (kappaHarm) * T_Bathx2 !+ boundray_term_vel(nx,iy,iz,T_Bathx2)
121-
end if
122-
end if
123-
end if
121+
end if
122+
end if
123+
end if clyBC
124124
end if
125125

126126
if (.not. Periodicy) then

‎src/heatflow/mod_evolve.f90‎

Lines changed: 110 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ module evolution
2727
use solver, only: linbcg
2828
use globe_data, only: Temp_p, Temp_pp, inverse_time, heat, lin_rhoc
2929
use globe_data, only: acsr, ja, ia
30+
use globe_data, only: heated_volume
3031
use heating, only: heater
3132
use boundary_vector, only: boundary
3233
use cattaneo, only: S_catS
@@ -117,6 +118,67 @@ subroutine simulate(itime)
117118
!---------------------------------------------
118119
! Construct S vector
119120
!---------------------------------------------
121+
122+
! COMPREHENSIVE DEBUG: Dump all quantities for radial cross-section at iy=16
123+
if (itime .le. 2) then
124+
block
125+
integer(int12) :: dbg_ix, dbg_idx, dbg_k
126+
real(real12) :: dbg_Ax, dbg_rowsum
127+
write(*,*) ''
128+
write(*,*) '=========================================================='
129+
write(*,'(A,I6)') ' DEBUG TIMESTEP ', itime
130+
write(*,*) '=========================================================='
131+
write(*,*) 'inverse_time =', inverse_time
132+
write(*,*) 'heated_volume =', heated_volume
133+
write(*,*) 'sum(Q) =', sum(Q), ' sum(Qdens) =', sum(Qdens)
134+
write(*,*) 'count(Qdens/=0) =', count(Qdens .ne. 0.0_real12)
135+
write(*,*) ''
136+
write(*,*) '--- Radial cross-section at iy=16, iz=1 ---'
137+
write(*,'(A6,A12,A12,A14,A14,A14,A14,A14)') &
138+
'ix', 'mat_id', 'kappa', 'rhoCp', 'Temp_p', 'B(I)', 'Qdens(I)', 'S(I)'
139+
do dbg_ix = 1, nx
140+
dbg_idx = dbg_ix + (16-1)*nx ! 1D index for (ix, iy=16, iz=1)
141+
write(*,'(I6,I12,ES12.4,ES14.6,ES14.6,ES14.6,ES14.6,ES14.6)') &
142+
dbg_ix, grid(dbg_ix,16,1)%imaterial_type, &
143+
grid(dbg_ix,16,1)%kappa, &
144+
lin_rhoc(dbg_idx), &
145+
Temp_p(dbg_idx), &
146+
B(dbg_idx), Qdens(dbg_idx), S(dbg_idx)
147+
end do
148+
write(*,*) ''
149+
write(*,*) '--- H-matrix rows for iy=16 (radial): row_sum and entries ---'
150+
do dbg_ix = 1, nx
151+
dbg_idx = dbg_ix + (16-1)*nx
152+
dbg_rowsum = 0.0_real12
153+
do dbg_k = ia(dbg_idx), ia(dbg_idx+1)-1
154+
dbg_rowsum = dbg_rowsum + acsr(dbg_k)
155+
end do
156+
! Compute A*Temp_p for this row (matrix-vector product)
157+
dbg_Ax = 0.0_real12
158+
do dbg_k = ia(dbg_idx), ia(dbg_idx+1)-1
159+
dbg_Ax = dbg_Ax + acsr(dbg_k) * Temp_p(ja(dbg_k))
160+
end do
161+
write(*,'(A,I3,A,ES14.6,A,ES14.6,A,ES14.6)') &
162+
' ix=', dbg_ix, &
163+
' row_sum=', dbg_rowsum, &
164+
' H*Tp=', dbg_Ax, &
165+
' S=', S(dbg_idx)
166+
end do
167+
write(*,*) ''
168+
write(*,*) '--- Heater region at iy=32, iz=1 ---'
169+
write(*,'(A6,A12,A14,A14,A14,A14)') &
170+
'ix', 'iheater', 'Temp_p', 'B(I)', 'Qdens(I)', 'S(I)'
171+
do dbg_ix = 1, min(10, nx)
172+
dbg_idx = dbg_ix + (32-1)*nx
173+
write(*,'(I6,I12,ES14.6,ES14.6,ES14.6,ES14.6)') &
174+
dbg_ix, grid(dbg_ix,32,1)%iheater, &
175+
Temp_p(dbg_idx), B(dbg_idx), Qdens(dbg_idx), S(dbg_idx)
176+
end do
177+
write(*,*) '=========================================================='
178+
write(*,*) ''
179+
end block
180+
end if
181+
120182
if ( iSteady .eq. 0 ) then
121183
S = - inverse_time * Temp_p * lin_rhoc - Qdens - B
122184
if (IVERB .gt. 3) then
@@ -192,12 +254,38 @@ subroutine simulate(itime)
192254

193255
call solve_petsc_csr(NA32, ia32, ja32, acsr, S, x, tol, itmax)
194256

195-
! Debug: Print solution statistics and verify solution
196-
if (IVERB .gt. 3) then
197-
write(*,*) "Solution x after PETSc: min=", minval(x), " max=", maxval(x), " avg=", sum(x)/size(x)
198-
write(*,*) "Temperature change: avg(x-Temp_p)=", sum(x-Temp_p)/size(x)
199-
write(*,*) "Max temperature change: ", maxval(abs(x-Temp_p))
200-
write(*,*) "=============================================="
257+
! POST-SOLVE DEBUG: Show solution and residual for radial cross-section
258+
if (itime .le. 2) then
259+
block
260+
integer(int12) :: dbg_ix, dbg_idx, dbg_k
261+
real(real12) :: dbg_Ax, dbg_resid
262+
write(*,*) ''
263+
write(*,*) '--- POST-SOLVE: Solution at iy=16, iz=1 ---'
264+
write(*,'(A6,A14,A14,A14,A14)') &
265+
'ix', 'Temp_p(old)', 'x(new)', 'deltaT', 'residual'
266+
do dbg_ix = 1, nx
267+
dbg_idx = dbg_ix + (16-1)*nx
268+
! Compute H*x for this row (should equal S)
269+
dbg_Ax = 0.0_real12
270+
do dbg_k = ia(dbg_idx), ia(dbg_idx+1)-1
271+
dbg_Ax = dbg_Ax + acsr(dbg_k) * x(ja(dbg_k))
272+
end do
273+
dbg_resid = dbg_Ax - S(dbg_idx)
274+
write(*,'(I6,ES14.6,ES14.6,ES14.6,ES14.6)') &
275+
dbg_ix, Temp_p(dbg_idx), x(dbg_idx), &
276+
x(dbg_idx) - Temp_p(dbg_idx), dbg_resid
277+
end do
278+
write(*,*) ''
279+
write(*,*) '--- POST-SOLVE: Heater region iy=32 ---'
280+
write(*,'(A6,A14,A14,A14)') 'ix', 'Temp_p(old)', 'x(new)', 'deltaT'
281+
do dbg_ix = 1, min(10, nx)
282+
dbg_idx = dbg_ix + (32-1)*nx
283+
write(*,'(I6,ES14.6,ES14.6,ES14.6)') &
284+
dbg_ix, Temp_p(dbg_idx), x(dbg_idx), &
285+
x(dbg_idx) - Temp_p(dbg_idx)
286+
end do
287+
write(*,*) '=========================================================='
288+
end block
201289
end if
202290

203291
! Note: Don't deallocate ia32, ja32 - keep them for next time step
@@ -227,6 +315,22 @@ subroutine simulate(itime)
227315
Temp_pp = Temp_p
228316
Temp_p = x
229317

318+
! DEBUG: Verify Temp_p after assignment
319+
if (itime .le. 2) then
320+
block
321+
integer(int12) :: dbg_ix2, dbg_idx2
322+
write(*,*) ''
323+
write(*,'(A,I6)') ' === VERIFY Temp_p AFTER ASSIGNMENT, itime=', itime
324+
write(*,'(A6,A14,A14)') 'ix', 'Temp_p(1D)', 'x(1D)'
325+
do dbg_ix2 = 1, nx
326+
dbg_idx2 = dbg_ix2 + (16-1)*nx
327+
write(*,'(I6,ES14.6,ES14.6)') &
328+
dbg_ix2, Temp_p(dbg_idx2), x(dbg_idx2)
329+
end do
330+
write(*,*) '=== END VERIFY ==='
331+
end block
332+
end if
333+
230334
! if (TempDepProp .eq. 1) then
231335
! CALL ChangeProp()
232336
! end if

‎src/heatflow/mod_hmatrix.f90‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -156,6 +156,8 @@ function hmatrixfunc(i, j) result(H)
156156
r_outer = real(x,real12)
157157
A = A * r_inner / r_center
158158
B = B * r_outer / r_center
159+
F = F * 0.0_real12
160+
G = G * 0.0_real12
159161
end if
160162

161163
! Determine the value of H based on the relationship between i and j

‎src/heatflow/mod_inputs.f90‎

Lines changed: 76 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -393,6 +393,34 @@ subroutine check_param(readvar,n)
393393
end if
394394

395395

396+
!------------------------------------------------------------------------------------
397+
! Cylindrical grid: only kappaBoundy1, kappaBoundNy, kappaBoundNr are needed.
398+
! Map them to the Cartesian variables and mark the others as satisfied.
399+
!------------------------------------------------------------------------------------
400+
if (CylindricalGrid) then
401+
! ix=1 is reflection axis -> zero flux (already handled in mod_boundary)
402+
kappaBoundx1 = 0.0
403+
readvar(9) = 1
404+
! ix=nx outer radius -> uses kappaBoundNr directly in mod_boundary
405+
! but also set kappaBoundNx so downstream code doesn't complain
406+
kappaBoundNx = kappaBoundNr
407+
readvar(29) = 1
408+
! nz=1 in cylindrical mode -> z boundaries are unused
409+
kappaBoundz1 = 0.0
410+
kappaBoundNz = 0.0
411+
readvar(11) = 1
412+
readvar(31) = 1
413+
! T_Bath mapping: x boundaries use cylindrical values
414+
T_Bathx1 = T_Bath ! reflection, not used
415+
T_Bathx2 = T_BathNr
416+
T_Bathz1 = T_Bath ! not used
417+
T_Bathz2 = T_Bath ! not used
418+
readvar(20) = 1 ! T_Bathx1
419+
readvar(21) = 1 ! T_Bathx2
420+
readvar(24) = 1 ! T_Bathz1
421+
readvar(25) = 1 ! T_Bathz2
422+
end if
423+
396424
PB:if ((.not. Periodicx).or.(.not. Periodicy).or.(.not. Periodicz)) then
397425
!------------------------------------------------------------------------------------
398426
! Error about missing kappa bound
@@ -423,28 +451,51 @@ subroutine check_param(readvar,n)
423451
write(6,'(A)') ' --- Warning in subroutine "check_param" ---'
424452
write(6,'(A)') ' --- Warning: KappaBound is not set ---'
425453
readvar(28) = 1
426-
427454
end if WarKBO
428455
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
429456
!------------------------------------------------------------------------------------
457+
! warning about missing kappa bound in cylindrical case. reassine to Kappabound
458+
!------------------------------------------------------------------------------------
459+
if (CylindricalGrid) then
460+
WarcylKB:if ( (readvar(49).eq.0) .and. (readvar(28) .eq. 1) ) then
461+
write(6,*)
462+
write(6,'(A43)') '###############################'
463+
write(6,'(A43)') '########## Warning ##########'
464+
write(6,'(A43)') '###############################'
465+
write(6,*)
466+
write(6,'(A)') ' --- Warning in subroutine "check_param" ---'
467+
write(6,'(A)') ' --- Warning: KappaBound is not set ---'
468+
readvar(49) = 1
469+
kappaBoundNr = KappaBound
470+
kappaBoundy1 = KappaBound
471+
kappaBoundNy = KappaBound
472+
kappaBoundz1 = 0
473+
kappaBoundNz = 0
474+
kappaBoundx1 = 0
475+
kappaBoundNx = 0
476+
end if WarcylKB
477+
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
478+
!------------------------------------------------------------------------------------
430479
! warning about missing kappa bound. reassine to Kappabound
431480
!------------------------------------------------------------------------------------
432-
WarKB:if (((any(readvar(9:11).eq.0)) .or. any(readvar(29:31).eq.0)) &
433-
.and. (readvar(28) .eq. 1) )then
434-
write(6,*)
435-
write(6,'(A43)') '###############################'
436-
write(6,'(A43)') '########## WARNING ##########'
437-
write(6,'(A43)') '###############################'
438-
write(6,*)
439-
write(6,'(A)') ' --- Warning in subroutine "check_param" ---'
440-
write(6,'(A)') ' --- WARNING: KappaBoundx,y,z not set, using KappaBound ---'
441-
kappaBoundx1 = KappaBound
442-
kappaBoundy1 = KappaBound
443-
kappaBoundz1 = KappaBound
444-
kappaBoundNx = KappaBound
445-
kappaBoundNy = KappaBound
446-
kappaBoundNz = KappaBound
447-
end if WarKB
481+
else
482+
WarKB:if (((any(readvar(9:11).eq.0)) .or. any(readvar(29:31).eq.0)) &
483+
.and. (readvar(28) .eq. 1) )then
484+
write(6,*)
485+
write(6,'(A43)') '###############################'
486+
write(6,'(A43)') '########## WARNING ##########'
487+
write(6,'(A43)') '###############################'
488+
write(6,*)
489+
write(6,'(A)') ' --- Warning in subroutine "check_param" ---'
490+
write(6,'(A)') ' --- WARNING: KappaBoundx,y,z not set, using KappaBound ---'
491+
kappaBoundx1 = KappaBound
492+
kappaBoundy1 = KappaBound
493+
kappaBoundz1 = KappaBound
494+
kappaBoundNx = KappaBound
495+
kappaBoundNy = KappaBound
496+
kappaBoundNz = KappaBound
497+
end if WarKB
498+
end if
448499
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
449500
elseif (((any(readvar(9:11).eq.0)) .or. any(readvar(29:31).eq.0)) &
450501
.or. (readvar(28) .gt. 0) ) then
@@ -655,7 +706,12 @@ subroutine read_system(unit)
655706
allocate(temp(nx))
656707
! read mesh volume dimessions
657708
read(unit,'(A)',iostat=Reason) buffer
658-
read(buffer,*) Lx, Ly, Lz
709+
if (CylindricalGrid) then
710+
read(buffer,*) Lx, Ly
711+
Lz = 1.0 !dummy value
712+
else
713+
read(buffer,*) Lx, Ly, Lz
714+
end if
659715
grid(:,:,:)%Length(1)=Lx/real(nx)
660716
grid(:,:,:)%Length(2)=Ly/real(ny)
661717
grid(:,:,:)%Length(3)=Lz/real(nz)
@@ -676,8 +732,8 @@ subroutine read_system(unit)
676732
r_in = real(ix - 1) * dr
677733
r_out = real(ix) * dr
678734
r_mid = 0.5_real12 * (r_in + r_out)
679-
grid(ix,:,:)%volume = pi * (r_out**2 - r_in**2) * grid(ix,1,1)%Length(2) &
680-
* grid(ix,1,1)%Length(3)
735+
grid(ix,:,:)%volume = pi * (r_out**2 - r_in**2) &
736+
* grid(ix,1,1)%Length(2)
681737
end do
682738
write(6,*) 'Cylindrical grid enabled: x=radial, y=axial, nz=1'
683739
end if

0 commit comments

Comments
 (0)