Skip to content

Commit 3b4da8a

Browse files
HarryMcleanHarryMclean
committed
Working Bicgstab
1 parent 5a375f5 commit 3b4da8a

6 files changed

Lines changed: 81 additions & 21 deletions

File tree

‎Makefile‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -42,9 +42,11 @@ MKL_INCLUDE_DIR = $(MKLROOT)/include
4242
#MKL_FLAGS = -L$(MKL_LIB_DIR) -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
4343
MKL_FLAGS = -L$(MKL_LIB_DIR) -lmkl_gf_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl
4444

45-
FFLAGS = -O3 -I$(MKL_INCLUDE_DIR)
45+
#FFLAGS = -O3 -I$(MKL_INCLUDE_DIR)
4646
MODULEFLAGS = -J$(BUILD_DIR)
4747
FC = gfortran
48+
NCORES := $(shell nproc)
49+
FFLAGS = -O3 -fopenmp -I$(MKL_INCLUDE_DIR)
4850

4951
##########################################
5052
# TARGETS
@@ -72,6 +74,10 @@ $(BUILD_DIR)/%.o: $(SRC_DIR)/%.f90 | $(BUILD_DIR)
7274
$(programs): $(OBJS) | $(BIN_DIR)
7375
$(FC) -O3 -fopenmp $(OBJS) -o $@ $(MKL_FLAGS)
7476

77+
.PHONY: run
78+
run: $(programs)
79+
OMP_NUM_THREADS=$(NCORES) MKL_NUM_THREADS=$(NCORES) MKL_DYNAMIC=FALSE OMP_PROC_BIND=spread OMP_PLACES=cores ./bin/$(NAME)
80+
7581
debug: FFLAGS = -O0 -Wall -g -ffpe-trap=invalid,zero,overflow,underflow -fbacktrace -fcheck=all -fbounds-check -I$(MKL_INCLUDE_DIR)
7682
debug: $(OBJS) | $(BIN_DIR)
7783
$(FC) $(FFLAGS) $(OBJS) -o $(programs) $(MKL_FLAGS)

‎src/heatflow/mod_boundary.f90‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -84,7 +84,7 @@ subroutine boundary(B)
8484
kappaHarm = (2*kappa*kappaBoundx1/(kappa+kappaBoundx1)) / &
8585
(grid(ix, iy, iz)%Length(1)**2)
8686
if (kappa .ne. kappaBoundx1) kappaHarm = kappaHarm*BR
87-
B(I) = B(I) + (kappaHarm) * T_Bathx1 + boundray_term_vel(1_int12,iy,iz,T_Bathx1)
87+
B(I) = B(I) + (kappaHarm) * T_Bathx1 !+ boundray_term_vel(1_int12,iy,iz,T_Bathx1)
8888
end if
8989
end if
9090
if (ix .eq. nx) then
@@ -94,7 +94,7 @@ subroutine boundary(B)
9494
kappaHarm = (2*kappa*kappaBoundNx/(kappa+kappaBoundNx)) / &
9595
(grid(ix, iy, iz)%Length(1)**2)
9696
if (kappa .ne. kappaBoundNx) kappaHarm = kappaHarm*BR
97-
B(I) = B(I) + (kappaHarm) * T_Bathx2 + boundray_term_vel(nx,iy,iz,T_Bathx2)
97+
B(I) = B(I) + (kappaHarm) * T_Bathx2 !+ boundray_term_vel(nx,iy,iz,T_Bathx2)
9898
end if
9999
end if
100100
end if
@@ -107,7 +107,7 @@ subroutine boundary(B)
107107
kappaHarm = (2*kappa*kappaBoundy1/(kappa+kappaBoundy1)) / &
108108
(grid(ix, iy, iz)%Length(2)**2)
109109
if (kappa .ne. kappaBoundy1) kappaHarm = kappaHarm*BR
110-
B(I) = B(I) + (kappaHarm) * T_Bathy1 + boundray_term_vel(ix,1_int12,iz,T_Bathy1)
110+
B(I) = B(I) + (kappaHarm) * T_Bathy1 !+ boundray_term_vel(ix,1_int12,iz,T_Bathy1)
111111
end if
112112
end if
113113
if (iy .eq. ny) then
@@ -117,7 +117,7 @@ subroutine boundary(B)
117117
kappaHarm = (2*kappa*kappaBoundNy/(kappa+kappaBoundNy)) / &
118118
(grid(ix, iy, iz)%Length(2)**2)
119119
if (kappa .ne. kappaBoundNy) kappaHarm = kappaHarm*BR
120-
B(I) = B(I) + (kappaHarm) * T_Bathy2 + boundray_term_vel(ix,ny,iz,T_Bathy2)
120+
B(I) = B(I) + (kappaHarm) * T_Bathy2 !+ boundray_term_vel(ix,ny,iz,T_Bathy2)
121121
end if
122122
end if
123123
end if
@@ -130,7 +130,7 @@ subroutine boundary(B)
130130
kappaHarm = (2*kappa*kappaBoundz1/(kappa+kappaBoundz1)) / &
131131
(grid(ix, iy, iz)%Length(3)**2)
132132
if (kappa .ne. kappaBoundz1) kappaHarm = kappaHarm*BR
133-
B(I) = B(I) + (kappaHarm) * T_Bathz1 + boundray_term_vel(ix,iy,1_int12,T_Bathz1)
133+
B(I) = B(I) + (kappaHarm) * T_Bathz1 !+ boundray_term_vel(ix,iy,1_int12,T_Bathz1)
134134
end if
135135
end if
136136
if (iz .eq. nz) then
@@ -140,7 +140,7 @@ subroutine boundary(B)
140140
kappaHarm = (2*kappa*kappaBoundNz/(kappa+kappaBoundNz)) / &
141141
(grid(ix, iy, iz)%Length(3)**2)
142142
if (kappa .ne. kappaBoundNz) kappaHarm = kappaHarm*BR
143-
B(I) = B(I) + (kappaHarm) * T_Bathz2 + boundray_term_vel(ix,iy,nz,T_Bathz2)
143+
B(I) = B(I) + (kappaHarm) * T_Bathz2 !+ boundray_term_vel(ix,iy,nz,T_Bathz2)
144144
end if
145145
end if
146146
end if

‎src/heatflow/mod_evolve.f90‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -150,9 +150,9 @@ subroutine simulate(itime)
150150
err=E
151151

152152

153-
! CALL bicgstab(acsr, ia, ja, S, itmax, x, x0, iter)
153+
call bicgstab(acsr, ia, ja, S, itmax, Temp_p, x, iter)
154154

155-
CALL solve_pardiso(acsr, S, ia, ja, x)
155+
! CALL solve_pardiso(acsr, S, ia, ja, x)
156156
! CALL linbcg(S,x,itol=int(itol,I4B),tol=tol, itmax=int(itmax,I4B), iter=iter, &
157157
! err=E)
158158

‎src/heatflow/mod_global.f90‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -14,8 +14,8 @@ module globe_data
1414
TYPE(sprs2_dp) :: ra !Techniqually rH
1515
TYPE(diag_sprs_dp) :: da !Techniqually dH
1616
real(real12), dimension(:), allocatable :: acsr
17-
integer, dimension(:), allocatable :: ja
18-
integer, dimension(:), allocatable :: ia
17+
integer(kind=8), dimension(:), allocatable :: ja
18+
integer(kind=8), dimension(:), allocatable :: ia
1919
character(len=1024) :: logname ! Log file name
2020

2121
end module globe_data

‎src/heatflow/mod_hmatrix.f90‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -149,7 +149,7 @@ function hmatrixfunc(i, j) result(H)
149149
H=0.0_real12
150150
else
151151
H = A ! X left neighbor (left cell interaction)
152-
H = H + calculate_convective_conductivity(xm, y, z, x, y, z)
152+
!H = H + calculate_convective_conductivity(xm, y, z, x, y, z)
153153
end if
154154
end if
155155

@@ -158,7 +158,7 @@ function hmatrixfunc(i, j) result(H)
158158
H=0.0_real12
159159
else
160160
H = B ! X right neighbor (right cell interaction)
161-
H = H + calculate_convective_conductivity(xp, y, z, x, y, z)
161+
!H = H + calculate_convective_conductivity(xp, y, z, x, y, z)
162162
end if
163163
end if
164164

@@ -167,15 +167,15 @@ function hmatrixfunc(i, j) result(H)
167167
H=0.0_real12
168168
else
169169
H = D ! Y down neighbor (down cell interaction)
170-
H = H + calculate_convective_conductivity(x, ym, z, x, y, z)
170+
!H = H + calculate_convective_conductivity(x, ym, z, x, y, z)
171171
end if
172172
end if
173173
if ((i-j) .eq. -nx) then
174174
if (y.eq.ny) then
175175
H=0.0_real12
176176
else
177177
H = E ! Y up neighbor (up cell interaction)
178-
H = H + calculate_convective_conductivity(x, yp, z, x, y, z)
178+
!H = H + calculate_convective_conductivity(x, yp, z, x, y, z)
179179
end if
180180
end if
181181

@@ -185,7 +185,7 @@ function hmatrixfunc(i, j) result(H)
185185
else
186186
!write(*,*) 'F this is forward (in) z',F
187187
H = F ! Z in neighbor (forward cell interaction) !!!Frank had this as G during testing
188-
H = H + calculate_convective_conductivity(x, y, zm, x, y, z)
188+
!H = H + calculate_convective_conductivity(x, y, zm, x, y, z)
189189
end if
190190
end if
191191

@@ -195,7 +195,7 @@ function hmatrixfunc(i, j) result(H)
195195
else
196196
!write(*,*) 'G this is backward (out) z?',G
197197
H = G ! Z out neighbor (backward cell interaction) !!!Frank Had this as F during testing
198-
H = H + calculate_convective_conductivity(x, y, zp, x, y, z)
198+
!H = H + calculate_convective_conductivity(x, y, zp, x, y, z)
199199
end if
200200
end if
201201

‎src/heatflow/mod_sparse_solver.f90‎

Lines changed: 58 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -56,8 +56,8 @@ subroutine coo2csr( nrow, &
5656
integer (kind=8), dimension(nnz), intent(in) :: ir
5757
integer (kind=8), dimension(nnz), intent(in) :: jc
5858
real(8), dimension(nnz), intent(out) :: acsr
59-
integer, dimension(nnz), intent(out) :: ja
60-
integer, dimension(nrow+1), intent(out) :: ia
59+
integer(kind=8), dimension(nnz), intent(out) :: ja
60+
integer(kind=8), dimension(nrow+1), intent(out) :: ia
6161

6262
! Local variables.
6363
integer (kind=8) :: i, iad, j, k, k0
@@ -240,7 +240,8 @@ subroutine solve_pardiso( acsr, &
240240
integer, dimension(:), allocatable :: iparm
241241
integer,dimension(1) :: idum
242242
real(8),dimension(1) :: ddum
243-
243+
integer :: badcol, missing_diag, k
244+
logical :: found
244245
n = size(b,1)
245246
nnz = size(acsr,1)
246247
nrhs = 1
@@ -284,7 +285,60 @@ subroutine solve_pardiso( acsr, &
284285
do i=1,64
285286
pt(i)%dummy = 0
286287
end do
287-
288+
289+
!---------------- CSR integrity / diagnostic checks ----------------
290+
291+
badcol = 0
292+
293+
write(*,*) 'PARDISO debug:'
294+
write(*,*) ' n =', n
295+
write(*,*) ' nnz =', nnz
296+
write(*,*) ' ia(1) =', ia(1), ' ia(n+1)=', ia(n+1), ' ia(n+1)-1=', ia(n+1)-1
297+
298+
if (ia(1) /= 1) stop 'ERROR: ia(1) must be 1'
299+
if (ia(n+1)-1 /= nnz) stop 'ERROR: ia end mismatch'
300+
301+
do i=1,n
302+
if (ia(i) > ia(i+1)) then
303+
write(*,*) 'Row pointer decreases at row', i
304+
stop 'ERROR: ia not monotone'
305+
end if
306+
end do
307+
308+
do k=1,nnz
309+
if (ja(k) < 1 .or. ja(k) > n) then
310+
badcol = badcol + 1
311+
if (badcol <= 10) write(*,*) 'Bad column index k=',k,' ja=',ja(k)
312+
end if
313+
end do
314+
if (badcol > 0) then
315+
write(*,*) 'Total bad columns =', badcol
316+
stop 'ERROR: invalid ja entries'
317+
end if
318+
319+
! Check each row has a diagonal and (optionally) detect duplicates
320+
missing_diag = 0
321+
do i=1,n
322+
found = .false.
323+
if (ia(i) < ia(i+1)) then
324+
! simple duplicate check (requires row segment unsorted ascending to be meaningful)
325+
do k = ia(i), ia(i+1)-1
326+
if (ja(k) == i) then
327+
if (acsr(k) == 0.0d0) then
328+
write(*,*) 'Zero diagonal at row', i
329+
stop 'ERROR: zero diagonal'
330+
end if
331+
found = .true.
332+
end if
333+
end do
334+
end if
335+
if (.not. found) then
336+
missing_diag = missing_diag + 1
337+
if (missing_diag <= 10) write(*,*) 'Missing diagonal at row', i
338+
end if
339+
end do
340+
if (missing_diag > 0) stop 'ERROR: missing diagonals'
341+
!-------------------------------------------------------------------
288342
phase = 11 ! Only reordering and symbolic factorization
289343

290344
call pardiso (pt,maxfct,mnum,mtype,phase,n,acsr,ia,ja, &

0 commit comments

Comments
 (0)