Skip to content

Commit 2bd3703

Browse files
HarryMcleanHarryMclean
committed
Switch to petsc
1 parent 3b4da8a commit 2bd3703

5 files changed

Lines changed: 235 additions & 79 deletions

File tree

‎Makefile‎

Lines changed: 116 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -1,89 +1,133 @@
11
####################################################################
2-
# 11 Jun 2024 #
2+
# HeatFlow Build System (MKL + optional PETSc + OpenMP)
33
####################################################################
4-
SHELL = /bin/sh
5-
PLAT = _linux
6-
7-
##########################################
8-
# CODE DIRECTORIES AND FILES
9-
##########################################
10-
mkfile_path := $(abspath $(firstword $(MAKEFILE_LIST)))
11-
mkfile_dir := $(dir $(mkfile_path))
12-
BIN_DIR := ./bin
13-
SRC_DIR := ./src
14-
BUILD_DIR = ./obj
15-
16-
SRCS := heatflow/mod_constants.f90 \
17-
heatflow/mod_constructions.f90 \
18-
heatflow/mod_SPtype.f90 \
19-
heatflow/mod_global.f90 \
20-
heatflow/mod_Sparse.f90 \
21-
heatflow/mod_inputs.f90 \
22-
heatflow/mod_material.f90 \
23-
heatflow/mod_hmatrix.f90 \
24-
heatflow/mod_init_evolve.f90 \
25-
heatflow/mkl_pardiso.f90 \
26-
heatflow/mod_sparse_solver.f90 \
27-
heatflow/mod_setup.f90 \
28-
heatflow/mod_boundary.f90 \
29-
heatflow/mod_heating.f90 \
30-
heatflow/mod_cattaneo.f90 \
31-
heatflow/mod_tempdep.f90 \
32-
heatflow/mod_evolve.f90 \
33-
heatflow/mod_output.f90 \
34-
heatflow.f90
354

36-
OBJS := $(addprefix $(BUILD_DIR)/,$(notdir $(SRCS:.f90=.o)))
37-
38-
# MKL configuration
39-
MKLROOT ?= /opt/intel/oneapi/mkl/latest
40-
MKL_LIB_DIR = $(MKLROOT)/lib/intel64
41-
MKL_INCLUDE_DIR = $(MKLROOT)/include
42-
#MKL_FLAGS = -L$(MKL_LIB_DIR) -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
43-
MKL_FLAGS = -L$(MKL_LIB_DIR) -lmkl_gf_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl
5+
SHELL = /bin/sh
6+
7+
# Directories
8+
SRC_DIR := ./src
9+
BUILD_DIR := ./obj
10+
BIN_DIR := ./bin
11+
12+
# Compiler
13+
FC := gfortran
14+
15+
# Core count
16+
NCORES := $(shell nproc)
17+
18+
# MKL
19+
MKLROOT ?= /opt/intel/oneapi/mkl/latest
20+
MKL_LIB_DIR := $(MKLROOT)/lib/intel64
21+
MKL_INCLUDE := $(MKLROOT)/include
22+
MKL_FLAGS := -L$(MKL_LIB_DIR) -lmkl_gf_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl
23+
24+
# PETSc (manual fallback if petsc-config missing)
25+
PETSC_PREFIX := /usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-real
26+
PETSC_FINCLUDE := /usr/share/petsc/3.15/include
27+
PETSC_AINCLUDE := $(PETSC_PREFIX)/include
28+
PETSC_LIBDIR := $(PETSC_PREFIX)/lib
29+
30+
PETSC_CONFIG := $(shell command -v petsc-config 2>/dev/null)
31+
ifeq ($(PETSC_CONFIG),)
32+
PETSC_INC := -I$(PETSC_FINCLUDE) -I$(PETSC_AINCLUDE)
33+
PETSC_LIB := -L$(PETSC_LIBDIR) -lpetsc
34+
PETSC_NOTE := (PETSc manual paths)
35+
else
36+
PETSC_INC := $(shell petsc-config --cflags)
37+
PETSC_LIB := $(shell petsc-config --libs)
38+
PETSC_NOTE := (petsc-config)
39+
endif
40+
41+
# Flags
42+
OPTFLAGS := -O3
43+
OMPFLAGS := -fopenmp
44+
WARNFLAGS := -Wall
45+
MODDIR_FLAG := -J$(BUILD_DIR)
46+
47+
FFLAGS := -cpp $(OPTFLAGS) $(OMPFLAGS) $(WARNFLAGS) -I$(MKL_INCLUDE) $(PETSC_INC) $(MODDIR_FLAG)
48+
DEBUGFLAGS := -cpp -O0 -g -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow,underflow -fbounds-check -I$(MKL_INCLUDE) $(PETSC_INC) $(MODDIR_FLAG)
49+
50+
# Program
51+
NAME := ThermalFlow.x
52+
TARGET := $(BIN_DIR)/$(NAME)
53+
54+
# Sources (module order)
55+
SRCS := \
56+
heatflow/mod_constants.f90 \
57+
heatflow/mod_constructions.f90 \
58+
heatflow/mod_SPtype.f90 \
59+
heatflow/mod_global.f90 \
60+
heatflow/mod_Sparse.f90 \
61+
heatflow/mod_inputs.f90 \
62+
heatflow/mod_material.f90 \
63+
heatflow/mod_hmatrix.f90 \
64+
heatflow/mod_init_evolve.f90 \
65+
heatflow/mkl_pardiso.f90 \
66+
heatflow/mod_sparse_solver.f90 \
67+
heatflow/mod_petsc_solver.f90 \
68+
heatflow/mod_boundary.f90 \
69+
heatflow/mod_heating.f90 \
70+
heatflow/mod_cattaneo.f90 \
71+
heatflow/mod_tempdep.f90 \
72+
heatflow/mod_evolve.f90 \
73+
heatflow/mod_output.f90 \
74+
heatflow/mod_setup.f90 \
75+
heatflow.f90
4476

45-
#FFLAGS = -O3 -I$(MKL_INCLUDE_DIR)
46-
MODULEFLAGS = -J$(BUILD_DIR)
47-
FC = gfortran
48-
NCORES := $(shell nproc)
49-
FFLAGS = -O3 -fopenmp -I$(MKL_INCLUDE_DIR)
50-
51-
##########################################
52-
# TARGETS
53-
##########################################
54-
NAME = ThermalFlow.x
55-
programs = $(BIN_DIR)/$(NAME)
77+
OBJS := $(addprefix $(BUILD_DIR)/,$(notdir $(SRCS:.f90=.o)))
5678

57-
.PHONY: all debug clean OMP
79+
.PHONY: all debug clean distclean run help show
5880

59-
all: $(programs)
81+
all: show $(TARGET)
6082

61-
$(BIN_DIR):
62-
mkdir -p $@
83+
show:
84+
@printf 'Building %s %s\n' '$(NAME)' '$(PETSC_NOTE)'
6385

64-
$(BUILD_DIR):
86+
$(BIN_DIR) $(BUILD_DIR):
6587
mkdir -p $@
6688

67-
# Pattern rule for compiling Fortran files
89+
# Compile module sources
6890
$(BUILD_DIR)/%.o: $(SRC_DIR)/heatflow/%.f90 | $(BUILD_DIR)
69-
$(FC) $(FFLAGS) $(MODULEFLAGS) -c $< -o $@
70-
71-
$(BUILD_DIR)/%.o: $(SRC_DIR)/%.f90 | $(BUILD_DIR)
72-
$(FC) $(FFLAGS) $(MODULEFLAGS) -c $< -o $@
91+
$(FC) $(FFLAGS) -c $< -o $@
7392

74-
$(programs): $(OBJS) | $(BIN_DIR)
75-
$(FC) -O3 -fopenmp $(OBJS) -o $@ $(MKL_FLAGS)
93+
# Main program
94+
$(BUILD_DIR)/heatflow.o: $(SRC_DIR)/heatflow.f90 | $(BUILD_DIR)
95+
$(FC) $(FFLAGS) -c $< -o $@
7696

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)
97+
# Link (single definition)
98+
$(TARGET): $(BIN_DIR) $(OBJS)
99+
$(FC) $(OPTFLAGS) $(OMPFLAGS) $(OBJS) -o $@ $(MKL_FLAGS) $(PETSC_LIB) -Wl,-rpath,$(PETSC_LIBDIR)
80100

81-
debug: FFLAGS = -O0 -Wall -g -ffpe-trap=invalid,zero,overflow,underflow -fbacktrace -fcheck=all -fbounds-check -I$(MKL_INCLUDE_DIR)
82-
debug: $(OBJS) | $(BIN_DIR)
83-
$(FC) $(FFLAGS) $(OBJS) -o $(programs) $(MKL_FLAGS)
101+
debug: FFLAGS = $(DEBUGFLAGS)
102+
debug: clean show $(TARGET)
84103

85-
OMP: $(programs)
86-
./util/DShell/omp_exec.sh
104+
run: $(TARGET)
105+
OMP_NUM_THREADS=$(NCORES) \
106+
MKL_NUM_THREADS=$(NCORES) \
107+
MKL_DYNAMIC=FALSE \
108+
OMP_PROC_BIND=spread \
109+
OMP_PLACES=cores \
110+
$< $(RUN_ARGS)
87111

88112
clean:
89-
rm -f $(BUILD_DIR)/*.o $(BUILD_DIR)/*.mod $(programs)
113+
@echo "[CLEAN] objects and modules"
114+
@rm -f $(BUILD_DIR)/*.o $(BUILD_DIR)/*.mod
115+
116+
distclean: clean
117+
@echo "[CLEAN] executable"
118+
@rm -f $(TARGET)
119+
120+
help:
121+
@echo "Targets:"
122+
@echo " make / make all - build optimized"
123+
@echo " make debug - debug build"
124+
@echo " make run - run with all cores"
125+
@echo " make clean - remove objects/modules"
126+
@echo " make distclean - remove executable"
127+
@echo "Variables:"
128+
@echo " RUN_ARGS='-ksp_type cg -pc_type gamg -ksp_rtol 1e-8 -ksp_monitor'"
129+
@echo "Parallel build: make -j$(NCORES)"
130+
131+
####################################################################
132+
# End
133+
####################################################################

‎src/heatflow.f90‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@ program HEATFLOW_V0_3
2828
use evolution, only: simulate
2929
use setup, only: set_global_variables
3030
use INITIAL, only: initial_evolve
31+
use petsc_solver, only: petsc_init, petsc_finalize
3132

3233
implicit none
3334
real(real12) :: cpustart, cpuend, cpustart2, progress
@@ -79,8 +80,11 @@ program HEATFLOW_V0_3
7980
! CALL initial_evolve to set systems initial Temperature conditions
8081
if (itime .eq. 1) CALL initial_evolve
8182

82-
! run the time evolution
83+
! run the time evolution
84+
CALL petsc_init()
8385
CALL simulate(itime)
86+
CALL petsc_finalize()
87+
8488

8589

8690
! Write results

‎src/heatflow/mod_evolve.f90‎

Lines changed: 18 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,8 @@ module evolution
3131
use boundary_vector, only: boundary
3232
use cattaneo, only: S_catS
3333
! use tempdep, only: ChangeProp
34-
use sparse_solver, only: bicgstab, solve_pardiso
34+
use petsc_solver, only: solve_petsc_csr
35+
3536
implicit none
3637

3738
private
@@ -51,7 +52,9 @@ subroutine simulate(itime)
5152
integer:: ncg, itol, itmax !, iss
5253
integer :: iter
5354
real(real12) :: e, err, tol
54-
55+
integer, allocatable :: ia32(:), ja32(:) ! 32-bit copies for PETSc
56+
integer :: NA32
57+
5558
!----------------------
5659
! Initialize vectors
5760
!----------------------
@@ -149,8 +152,19 @@ subroutine simulate(itime)
149152
iter= 0
150153
err=E
151154

152-
153-
call bicgstab(acsr, ia, ja, S, itmax, Temp_p, x, iter)
155+
print *, "Calling solver..."
156+
! call bicgstab(acsr, ia, ja, S, itmax, Temp_p, x, iter)
157+
if (.not. allocated(ia32)) then
158+
allocate(ia32(size(ia)), ja32(size(ja)))
159+
ia32 = int(ia, kind=kind(ia32))
160+
ja32 = int(ja, kind=kind(ja32))
161+
end if
162+
NA32 = int(NA, kind=kind(NA32))
163+
allocate(x(NA))
164+
call solve_petsc_csr(NA32, ia32, ja32, acsr, S, x, tol, itmax)
165+
if (allocated(ia32)) deallocate(ia32, ja32)
166+
167+
print *, "Solver finished."
154168

155169
! CALL solve_pardiso(acsr, S, ia, ja, x)
156170
! CALL linbcg(S,x,itol=int(itol,I4B),tol=tol, itmax=int(itmax,I4B), iter=iter, &

‎src/heatflow/mod_petsc_solver.f90‎

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
module petsc_solver
2+
#include "petsc/finclude/petscsys.h"
3+
#include "petsc/finclude/petscksp.h"
4+
use petscksp
5+
implicit none
6+
private
7+
public :: petsc_init, petsc_finalize, solve_petsc_csr
8+
9+
contains
10+
11+
subroutine petsc_init()
12+
integer :: ierr
13+
call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
14+
end subroutine petsc_init
15+
16+
subroutine petsc_finalize()
17+
integer :: ierr
18+
call PetscFinalize(ierr)
19+
end subroutine petsc_finalize
20+
21+
subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
22+
integer, intent(in) :: n
23+
integer, intent(in) :: ia(:), ja(:)
24+
real(8), intent(in) :: aval(:), b(:)
25+
real(8), intent(inout) :: x(:)
26+
real(8), intent(in) :: rtol
27+
integer, intent(in) :: maxit
28+
Mat :: A
29+
Vec :: bb, xx
30+
KSP :: ksp
31+
PC :: pc
32+
integer :: ierr, i, row_nz, start_k
33+
integer, allocatable :: cols0(:), idx(:)
34+
real(8), allocatable :: vals(:)
35+
real(8), pointer :: xptr(:)
36+
37+
if (size(ia) /= n+1) stop 'solve_petsc_csr: ia size mismatch'
38+
if (size(b) /= n .or. size(x) /= n) stop 'solve_petsc_csr: vector size mismatch'
39+
40+
! Create matrix with an estimated 7 nonzeros/row (adjust if needed)
41+
call MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 7, PETSC_NULL_INTEGER, A, ierr)
42+
43+
do i = 1, n
44+
row_nz = ia(i+1) - ia(i)
45+
if (row_nz > 0) then
46+
start_k = ia(i)
47+
allocate(cols0(row_nz), vals(row_nz))
48+
cols0 = ja(start_k:start_k+row_nz-1) - 1 ! zero-based
49+
vals = aval(start_k:start_k+row_nz-1)
50+
call MatSetValues(A, 1, (/i-1/), row_nz, cols0, vals, INSERT_VALUES, ierr)
51+
deallocate(cols0, vals)
52+
end if
53+
end do
54+
call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)
55+
call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)
56+
57+
! Create vectors
58+
call VecCreateSeq(PETSC_COMM_SELF, n, bb, ierr)
59+
call VecCreateSeq(PETSC_COMM_SELF, n, xx, ierr)
60+
61+
! Set RHS and initial guess
62+
allocate(idx(n))
63+
idx = [(i-1, i=1,n)]
64+
call VecSetValues(bb, n, idx, b, INSERT_VALUES, ierr)
65+
call VecAssemblyBegin(bb,ierr); call VecAssemblyEnd(bb,ierr)
66+
67+
call VecSetValues(xx, n, idx, x, INSERT_VALUES, ierr)
68+
call VecAssemblyBegin(xx,ierr); call VecAssemblyEnd(xx,ierr)
69+
70+
! KSP setup
71+
call KSPCreate(PETSC_COMM_SELF, ksp, ierr)
72+
call KSPSetOperators(ksp, A, A, ierr) ! 3-arg form (reuse automatically)
73+
call KSPGetPC(ksp, pc, ierr)
74+
call PCSetType(pc, PCILU, ierr) ! Override at runtime with -pc_type
75+
call KSPSetType(ksp, KSPCG, ierr) ! Use -ksp_type bcgs if not SPD
76+
call KSPSetTolerances(ksp, rtol, PETSC_DEFAULT_REAL, PETSC_DEFAULT_REAL, maxit, ierr)
77+
call KSPSetFromOptions(ksp, ierr)
78+
79+
call KSPSolve(ksp, bb, xx, ierr)
80+
81+
! Extract solution
82+
call VecGetArrayF90(xx, xptr, ierr)
83+
x = xptr
84+
call VecRestoreArrayF90(xx, xptr, ierr)
85+
86+
! Cleanup
87+
deallocate(idx)
88+
call KSPDestroy(ksp, ierr)
89+
call VecDestroy(bb, ierr)
90+
call VecDestroy(xx, ierr)
91+
call MatDestroy(A, ierr)
92+
end subroutine solve_petsc_csr
93+
94+
end module petsc_solver

‎src/heatflow/mod_sparse_solver.f90‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -139,7 +139,6 @@ subroutine bicgstab( acsr, &
139139
real(8), dimension(:), allocatable :: r, p, s, rst, temp1, temp2
140140

141141
n = size(b,1)
142-
143142
allocate(x(n))
144143
allocate(r(n))
145144
allocate(p(n))
@@ -150,6 +149,7 @@ subroutine bicgstab( acsr, &
150149

151150
call mkl_dcsrgemv("N",n,acsr,ia,ja,x,temp1)
152151

152+
print *, "Initial residual norm: ", norm2(b - temp1)
153153
r = b - temp1
154154

155155
call random_number(rst)
@@ -163,7 +163,7 @@ subroutine bicgstab( acsr, &
163163
delta0 = delta
164164

165165
do i = 1, maxiter
166-
166+
print *, "Iteration ", i
167167
if ( norm2(r) /= norm2(r) ) then
168168
write(*,'(a)') "Error in solver: residual NaN"
169169
exit

0 commit comments

Comments
 (0)