diff --git a/CHANGELOG b/CHANGELOG index 794b181..35a47e5 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -12,6 +12,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed ### Security +##[2.0.0] - 2026-03-26 +## Added +-Add Petsc solver +-Add new constant gradient boundary conditions +-Add new compressed output flag + + + + ##[1.0.0] - 2025-07-28 ## Added -Additional boundary conditions diff --git a/Makefile b/Makefile index 57e89be..d7715eb 100644 --- a/Makefile +++ b/Makefile @@ -1,70 +1,172 @@ #################################################################### -# 11 Jun 2024 # +# HeatFlow Build System (MKL + optional PETSc + OpenMP) #################################################################### -# -SHELL = /bin/sh -# -# The machine (platform) identifier to append to the library names -# -PLAT = _linux -# -# - - -########################################## -# CODE DIRECTORIES AND FILES -########################################## -mkfile_path := $(abspath $(firstword $(MAKEFILE_LIST))) -mkfile_dir := $(dir $(mkfile_path)) -BIN_DIR := ./bin -SRC_DIR := ./src -BUILD_DIR = ./obj - -SRCS := /heatflow/mod_constants.f90 \ - /heatflow/mod_constructions.f90 \ - /heatflow/mod_SPtype.f90 \ - /heatflow/mod_global.f90 \ - /heatflow/mod_Sparse.f90 \ - /heatflow/mod_inputs.f90 \ - /heatflow/mod_material.f90 \ - /heatflow/mod_hmatrix.f90 \ - /heatflow/mod_init_evolve.f90 \ - /heatflow/mod_setup.f90 \ - /heatflow/mod_boundary.f90 \ - /heatflow/mod_heating.f90 \ - /heatflow/mod_cattaneo.f90 \ - /heatflow/mod_tempdep.f90 \ - /heatflow/mod_evolve.f90 \ - /heatflow/mod_output.f90 \ - heatflow.f90 -OBJS := $(addprefix $(SRC_DIR)/,$(SRCS)) - - -FFLAGS = -O3 -MODULEFLAGS = -J -FC = gfortran - -########################################## -# LIBRARY SECTION -########################################## -MKLROOT?="/usr/local/intel/parallel_studio_xe_2017/compilers_and_libraries_2017/linux/mkl/lib/intel64_lin" - - -NAME = ThermalFlow.x -programs = $(BIN_DIR)/$(NAME) -all: $(programs) - -$(BIN_DIR): - mkdir -p $@ -$(BUILD_DIR): +SHELL = /bin/sh + +# Directories +SRC_DIR := ./src +BUILD_DIR := ./obj +BIN_DIR := ./bin + +# Compiler (prefer system MPI wrapper for PETSc builds, allow user override) +SYSTEM_PATH := PATH=/usr/bin:/bin +ifeq ($(origin FC), default) +ifneq ($(wildcard /usr/bin/mpifort),) +FC := env $(SYSTEM_PATH) /usr/bin/mpifort +else ifneq ($(wildcard /usr/bin/gfortran),) +FC := env $(SYSTEM_PATH) /usr/bin/gfortran +else +FC := gfortran +endif +endif + +# Core count +NCORES := $(shell nproc) + +# Detect conda environment for BLAS/LAPACK (fallback if no system libs) +CONDA_PREFIX ?= $(shell conda info --base 2>/dev/null || echo /home/hm556/miniforge3) + +# PETSc (discover dynamically when possible) +ifneq ($(wildcard /usr/bin/pkg-config),) +PKG_CONFIG := env $(SYSTEM_PATH) /usr/bin/pkg-config +else +PKG_CONFIG := pkg-config +endif +PETSC_PKG_CFLAGS := $(shell $(PKG_CONFIG) --cflags petsc 2>/dev/null || $(PKG_CONFIG) --cflags PETSc 2>/dev/null) +PETSC_PKG_LIBS := $(shell $(PKG_CONFIG) --libs petsc 2>/dev/null || $(PKG_CONFIG) --libs PETSc 2>/dev/null) + +ifdef PETSC_DIR +PETSC_DIR_INC := -I$(PETSC_DIR)/include +ifdef PETSC_ARCH +PETSC_DIR_INC += -I$(PETSC_DIR)/$(PETSC_ARCH)/include +PETSC_DIR_LIB := -L$(PETSC_DIR)/$(PETSC_ARCH)/lib -Wl,-rpath,$(PETSC_DIR)/$(PETSC_ARCH)/lib +endif +endif + +ifeq ($(strip $(PETSC_PKG_CFLAGS)),) +PETSC_INC := $(PETSC_DIR_INC) +PETSC_LIB := $(PETSC_DIR_LIB) -lpetsc +PETSC_NOTE := (PETSc from PETSC_DIR/PETSC_ARCH) +else +PETSC_INC := $(PETSC_PKG_CFLAGS) +PETSC_LIB := $(PETSC_PKG_LIBS) +PETSC_NOTE := (PETSc via pkg-config) +endif + +ifeq ($(strip $(PETSC_INC)),) +PETSC_INC := -I/usr/share/petsc/3.15/include -I/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-real/include +PETSC_LIB := -L/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-real/lib -lpetsc -Wl,-rpath,/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-real/lib +PETSC_NOTE := (legacy PETSc 3.15 fallback) +endif + +# BLAS/LAPACK backend +OPENBLAS_LIBS := $(shell $(PKG_CONFIG) --libs openblas 2>/dev/null) +CONDA_BLAS_LIB := $(firstword $(wildcard $(CONDA_PREFIX)/lib/libblas.so.3 $(CONDA_PREFIX)/lib/libblas.so)) +CONDA_LAPACK_LIB := $(firstword $(wildcard $(CONDA_PREFIX)/lib/liblapack.so.3 $(CONDA_PREFIX)/lib/liblapack.so)) +ifeq ($(strip $(OPENBLAS_LIBS)),) +ifneq ($(strip $(CONDA_BLAS_LIB)$(CONDA_LAPACK_LIB)),) +BLAS_FLAGS := -Wl,--disable-new-dtags -Wl,-rpath,$(CONDA_PREFIX)/lib -Wl,--no-as-needed $(CONDA_BLAS_LIB) $(CONDA_LAPACK_LIB) -Wl,--as-needed -lpthread -lm +BLAS_NOTE := (OpenBLAS from CONDA_PREFIX) +else +BLAS_FLAGS := -lblas -llapack -lpthread -lm +BLAS_NOTE := (system BLAS/LAPACK fallback) +endif +else +BLAS_FLAGS := $(OPENBLAS_LIBS) -lpthread -lm +BLAS_NOTE := (OpenBLAS via pkg-config) +endif + +# Flags +OPTFLAGS := -O3 +OMPFLAGS := -fopenmp +WARNFLAGS := -Wall +MODDIR_FLAG := -J$(BUILD_DIR) + +FFLAGS := -cpp $(OPTFLAGS) $(OMPFLAGS) $(WARNFLAGS) $(PETSC_INC) $(MODDIR_FLAG) +DEBUGFLAGS := -cpp -O0 -g -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow,underflow -fbounds-check $(PETSC_INC) $(MODDIR_FLAG) + +# Program +NAME := ThermalFlow.x +TARGET := $(BIN_DIR)/$(NAME) + +# Sources (module order) +SRCS := \ + heatflow/mod_constants.f90 \ + heatflow/mod_constructions.f90 \ + heatflow/mod_SPtype.f90 \ + heatflow/mod_global.f90 \ + heatflow/mod_Sparse.f90 \ + heatflow/mod_inputs.f90 \ + heatflow/mod_material.f90 \ + heatflow/mod_hmatrix.f90 \ + heatflow/mod_init_evolve.f90 \ + heatflow/mod_petsc_solver.f90 \ + heatflow/mod_boundary.f90 \ + heatflow/mod_heating.f90 \ + heatflow/mod_cattaneo.f90 \ + heatflow/mod_tempdep.f90 \ + heatflow/mod_evolve.f90 \ + heatflow/mod_output.f90 \ + heatflow/mod_setup.f90 \ + heatflow.f90 + +OBJS := $(addprefix $(BUILD_DIR)/,$(notdir $(SRCS:.f90=.o))) + +.NOTPARALLEL: + +.PHONY: all debug clean distclean run help show + +all: show $(TARGET) + +show: + @printf 'Building %s %s %s\n' '$(NAME)' '$(PETSC_NOTE)' '$(BLAS_NOTE)' + +$(BIN_DIR) $(BUILD_DIR): mkdir -p $@ -$(programs) : $(OBJS) | $(BIN_DIR) $(BUILD_DIR) - $(FC) -O3 -fopenmp $(MODULEFLAGS) $(BUILD_DIR) $(OBJS) -o $@ +# Compile module sources +$(BUILD_DIR)/%.o: $(SRC_DIR)/heatflow/%.f90 | $(BUILD_DIR) + $(FC) $(FFLAGS) -c $< -o $@ + +# Main program +$(BUILD_DIR)/heatflow.o: $(SRC_DIR)/heatflow.f90 | $(BUILD_DIR) + $(FC) $(FFLAGS) -c $< -o $@ -debug : $(OBJS) - $(FC) -O0 -Wall -g -ffpe-trap=invalid,zero,overflow,underflow -fbacktrace -fcheck=all -fbounds-check $(MODULEFLAGS) $(BUILD_DIR) $(OBJS) -o $(programs) +# Link (single definition) +$(TARGET): $(BIN_DIR) $(OBJS) + $(FC) $(OPTFLAGS) $(OMPFLAGS) $(OBJS) -o $@ $(BLAS_FLAGS) $(PETSC_LIB) -OMP: $(programs) - ./util/DShell/omp_exec.sh +debug: FFLAGS = $(DEBUGFLAGS) +debug: clean show $(TARGET) + +run: $(TARGET) + mpiexec -n $(NCORES) \ + OMP_NUM_THREADS=1 \ + OPENBLAS_NUM_THREADS=1 \ + OMP_PROC_BIND=spread \ + OMP_PLACES=cores \ + $< $(RUN_ARGS) + +clean: + @echo "[CLEAN] objects and modules" + @rm -f $(BUILD_DIR)/*.o $(BUILD_DIR)/*.mod + +distclean: clean + @echo "[CLEAN] executable" + @rm -f $(TARGET) + +help: + @echo "Targets:" + @echo " make / make all - build optimized" + @echo " make debug - debug build" + @echo " make run - run distributed across all cores via MPI" + @echo " make clean - remove objects/modules" + @echo " make distclean - remove executable" + @echo "Variables:" + @echo " RUN_ARGS='-ksp_type cg -pc_type gamg -ksp_rtol 1e-8 -ksp_monitor'" + @echo "Parallel build: make -j$(NCORES)" + +#################################################################### +# End +#################################################################### \ No newline at end of file diff --git a/docs/USER_MANUAL.md b/docs/USER_MANUAL.md new file mode 100644 index 0000000..c622bc9 --- /dev/null +++ b/docs/USER_MANUAL.md @@ -0,0 +1,140 @@ +# HeatFlow User Manual + +This manual provides a concise guide to configuring and running simulations using the **HeatFlow** software. The software simulates heat transport using finite difference methods, primarily focusing on the Cattaneo (hyperbolic heat equation) and Fourier models. + +## Input Files + +The simulation is controlled by three main input files located in the `inputs/` directory: +1. **`param.in`**: Simulation parameters (time steps, flags, boundary conditions). +2. **`mat.in`**: Material properties. +3. **`system.in`**: Geometry and grid definition. + +### 1. `param.in` (Simulation Parameters) + +This file uses a `KEYWORD = VALUE` format. Comments can be added using `!`. + +#### General Settings +| Keyword | Type | Default | Description | +| :--- | :--- | :--- | :--- | +| `_RunName` | String | `default` | Name of the simulation run. | +| `IVERB` | Integer | `1` | Verbosity level (higher = more output). | +| `ntime` | Integer | `10` | Total number of time steps. | +| `time_step` | Double | `1.0` | Time step size. | +| `freq` | Double | `1.0` | Frequency of the heater. | +| `icattaneo` | Integer | `1` | Switch for Cattaneo term (`1` = On, `0` = Off/Fourier). | +| `isteady` | Integer | `0` | Steady state switch (`1` = Steady state, `0` = Transient). | +| `heattime` | Integer | `0` | Number of steps for which heating is applied (case 2). | +| `TempDepProp`| Integer | `0` | Flag for temperature dependent properties. | + +#### Boundary & Conditions +| Keyword | Type | Default | Description | +| :--- | :--- | :--- | :--- | +| `iboundary` | Integer | `1` | Boundary condition type. | +| `Periodic` | String | `''` | Periodic boundaries. Contains 'x', 'y', or 'z' (e.g., `'xy'`). | +| `kappaBound` | Double | `0.0` | Global boundary thermal conductivity (sets all planes). | +| `kappaBoundx1`...`z2` | Double | `0.0` | Specific boundary conductivity (e.g., `kappaBoundx1` for x=1 plane). | +| `T_System` | Double | `300.0` | Initial system temperature. | +| `T_Bath` | Double | - | Global bath temperature (sets all boundaries boundaries). | +| `T_Bathx1`...`z2` | Double | `T_Bath` | Specific boundary temperatures. | +| `T_BathCG` | Double | `0.0` | Constant gradient bath temperature. | +| `CG_dir` | String | `' '` | Direction for constant gradient (e.g., `'+x'`, `'-y'`). | +| `T_BathCC` | Logical| `F` | Scale constant gradient with DeltaT. | +| `BR` | Double | `1.0` | Bath Ratio (scaling factor). | + +#### Power +| Keyword | Type | Default | Description | +| :--- | :--- | :--- | :--- | +| `power_in` | Double | `0.0` | Power input for the heater. | + +#### Flags (Logical) +All flags default to `.False.`. Set to `.True.` (or `T`) to enable. +- `_Check_Sparse_Full`: Check if simulation is sparse or full. +- `_Check_Stability`: Perform stability check. +- `_Check_Steady_State`: Check for steady state convergence. +- `_WriteToTxt`: Enable writing output to text files. +- `_Percentage_Completion`: Show progress % in output. +- `_Test_Run`: Flag for test runs. +- `_InputTempDis`: Load initial temperature distribution from file. +- `_FullRestart`: Perform a full restart. + +#### Output Control +Defines the region of the grid to write to output. +| Keyword | Type | Default | Description | +| :--- | :--- | :--- | :--- | +| `write_every` | Integer | `1` | Write output every N steps. | +| `start_ix`, `end_ix` | Integer | `1`..`Nx` | X-range for output. | +| `start_iy`, `end_iy` | Integer | `1`..`Ny` | Y-range for output. | +| `start_iz`, `end_iz` | Integer | `1`..`Nz` | Z-range for output. | + +--- + +### 2. `mat.in` (Material Properties) + +Defines the physical properties for each material index used in the system. The file ends with a line containing `0`. + +**Format:** +``` + +keyword = value +... +0 +``` + +| Keyword | Description | +| :--- | :--- | +| `heat_capacity` | Specific heat capacity. | +| `kappa` | Thermal conductivity. | +| `rho` | Density. | +| `tau` | Relaxation time (for Cattaneo). | +| `em` | Emissivity / Parameter (usage depends on physics context). | +| `vel` | Velocity vector (3 components, e.g., `vel = 1.0 0.0 0.0`). | + +**Example:** +``` +1 +heat_capacity = 4200 +kappa = 0.541 +rho = 997 +tau = 1e-12 +vel = 0.0 0.0 0.0 +0 +``` + +--- + +### 3. `system.in` (Geometry/Mesh) + +Defines the simulation grid and the material distribution. + +**Structure:** +1. **Grid Dimensions**: `nx ny nz` +2. **Physical Dimensions**: `Lx Ly Lz` +3. **Grid Data**: A list of `MaterialID:HeaterID` for every cell. + +The file is read in the order: Z-planes, then Y-rows, then X-columns. +Each line in the file (after header) corresponds to one row (X-direction). + +**Example:** +``` +10 10 1 +0.01 0.01 0.001 + +! Z=1, Y=1 Row +1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 +! Z=1, Y=2 Row +1:0 1:0 ... +``` +- `1:0` means Material ID 1, Heater ID 0 (no heater). +- `1:1` means Material ID 1, Heater ID 1 (active heater). + +## Execution + +Ensure the `inputs/` directory exists with the three required files. Run the executable from the directory containing `inputs/`. + +```bash +./ThermalFlow.x +``` +or via `fpm`: +```bash +fpm run --profile release +``` diff --git a/fpm.toml b/fpm.toml index a1f0b5d..fcc5f3f 100644 --- a/fpm.toml +++ b/fpm.toml @@ -9,12 +9,24 @@ description = "A Fortran executable for modelling Cattaneo heat flow" [preprocess.cpp] suffixes = ["F90", "f90"] +[dependencies] +mpi = "*" + +[features] +petsc.build.external-modules = ["petscksp"] +petsc.build.link = ["petsc"] +petsc.flags = "-I../petsc/include -I../petsc/arch-darwin-c-debug/include" +petsc.link-time-flags = "-L../petsc/arch-darwin-c-debug/lib -Wl,-rpath,../petsc/arch-darwin-c-debug/lib" + +[profiles] +default = ["petsc"] + [library] source-dir="src" [fortran] -implicit-typing = false -implicit-external = false +implicit-typing = true +implicit-external = true source-form = "free" [[executable]] diff --git a/src/.vscode/settings.json b/src/.vscode/settings.json new file mode 100644 index 0000000..a8c2003 --- /dev/null +++ b/src/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "python-envs.defaultEnvManager": "ms-python.python:conda", + "python-envs.defaultPackageManager": "ms-python.python:conda", + "python-envs.pythonProjects": [] +} \ No newline at end of file diff --git a/src/heatflow.f90 b/src/heatflow.f90 index 4b4df6e..cc16531 100644 --- a/src/heatflow.f90 +++ b/src/heatflow.f90 @@ -28,11 +28,18 @@ program HEATFLOW_V0_3 use evolution, only: simulate use setup, only: set_global_variables use INITIAL, only: initial_evolve + use petsc_solver, only: petsc_init, petsc_finalize, petsc_is_root implicit none real(real12) :: cpustart, cpuend, cpustart2, progress integer(int12) :: itime + !-------------------------------------------------------------! + ! Initialize PETSc FIRST (before any other operations) ! + !-------------------------------------------------------------! + CALL petsc_init() + !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^! + !-------------------------------------------------------------! ! calculate the time to run full simulation ! !-------------------------------------------------------------! @@ -40,7 +47,7 @@ program HEATFLOW_V0_3 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^! ! give feedback to user that code has begun - write(*,*) 'Setup initialising' + if (petsc_is_root()) write(*,*) 'Setup initialising' !-------------------------------------------------------------! ! Read parameters from input file and set global variables ...! @@ -59,14 +66,15 @@ program HEATFLOW_V0_3 ! give feedback to user that main simulation is begining - write(*,*) 'Setup complete, running simulation' + if (petsc_is_root()) write(*,*) 'Setup complete, running simulation' !-------------------------------------------------------------! ! run simulation for 'ntime' time steps ! !-------------------------------------------------------------! + do itime=1,ntime - if (iverb.eq.0) then + if (petsc_is_root() .and. iverb.eq.0) then if (Lpercentage) then progress = real(itime)/real(ntime)*100.0 write(*,'(A,A,F12.4,A)', advance = 'no') achar(13)& @@ -79,26 +87,29 @@ program HEATFLOW_V0_3 ! CALL initial_evolve to set systems initial Temperature conditions if (itime .eq. 1) CALL initial_evolve - ! run the time evolution + ! run the time evolution CALL simulate(itime) + ! Write results - CALL data_write(itime) - if (IVERB.ge.3) CALL final_print + if (petsc_is_root()) CALL data_write(itime) + if (petsc_is_root() .and. IVERB.ge.3) CALL final_print - end do + end do + CALL petsc_finalize() + !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^! !-------------------------------------------------------------! ! calculate end time and print to user ! !-------------------------------------------------------------! - CALL cpu_time(cpuend) - write(*,'(A,F12.6)') ' time=', cpuend-cpustart + CALL cpu_time(cpuend) + if (petsc_is_root()) write(*,'(A,F12.6)') ' time=', cpuend-cpustart !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^! ! give feedback to user that code has ended - write(*,*) 'all done' + if (petsc_is_root()) write(*,*) 'all done' end program HEATFLOW_V0_3 diff --git a/src/heatflow/mod_boundary.f90 b/src/heatflow/mod_boundary.f90 index d94a0bf..88ad7ba 100644 --- a/src/heatflow/mod_boundary.f90 +++ b/src/heatflow/mod_boundary.f90 @@ -84,7 +84,7 @@ subroutine boundary(B) kappaHarm = (2*kappa*kappaBoundx1/(kappa+kappaBoundx1)) / & (grid(ix, iy, iz)%Length(1)**2) if (kappa .ne. kappaBoundx1) kappaHarm = kappaHarm*BR - B(I) = B(I) + (kappaHarm) * T_Bathx1 + B(I) = B(I) + (kappaHarm) * T_Bathx1 !+ boundray_term_vel(1_int12,iy,iz,T_Bathx1) end if end if if (ix .eq. nx) then @@ -94,7 +94,7 @@ subroutine boundary(B) kappaHarm = (2*kappa*kappaBoundNx/(kappa+kappaBoundNx)) / & (grid(ix, iy, iz)%Length(1)**2) if (kappa .ne. kappaBoundNx) kappaHarm = kappaHarm*BR - B(I) = B(I) + (kappaHarm) * T_Bathx2 + B(I) = B(I) + (kappaHarm) * T_Bathx2 !+ boundray_term_vel(nx,iy,iz,T_Bathx2) end if end if end if @@ -107,7 +107,7 @@ subroutine boundary(B) kappaHarm = (2*kappa*kappaBoundy1/(kappa+kappaBoundy1)) / & (grid(ix, iy, iz)%Length(2)**2) if (kappa .ne. kappaBoundy1) kappaHarm = kappaHarm*BR - B(I) = B(I) + (kappaHarm) * T_Bathy1 + B(I) = B(I) + (kappaHarm) * T_Bathy1 !+ boundray_term_vel(ix,1_int12,iz,T_Bathy1) end if end if if (iy .eq. ny) then @@ -117,7 +117,7 @@ subroutine boundary(B) kappaHarm = (2*kappa*kappaBoundNy/(kappa+kappaBoundNy)) / & (grid(ix, iy, iz)%Length(2)**2) if (kappa .ne. kappaBoundNy) kappaHarm = kappaHarm*BR - B(I) = B(I) + (kappaHarm) * T_Bathy2 + B(I) = B(I) + (kappaHarm) * T_Bathy2 !+ boundray_term_vel(ix,ny,iz,T_Bathy2) end if end if end if @@ -130,7 +130,7 @@ subroutine boundary(B) kappaHarm = (2*kappa*kappaBoundz1/(kappa+kappaBoundz1)) / & (grid(ix, iy, iz)%Length(3)**2) if (kappa .ne. kappaBoundz1) kappaHarm = kappaHarm*BR - B(I) = B(I) + (kappaHarm) * T_Bathz1 + B(I) = B(I) + (kappaHarm) * T_Bathz1 !+ boundray_term_vel(ix,iy,1_int12,T_Bathz1) end if end if if (iz .eq. nz) then @@ -140,7 +140,7 @@ subroutine boundary(B) kappaHarm = (2*kappa*kappaBoundNz/(kappa+kappaBoundNz)) / & (grid(ix, iy, iz)%Length(3)**2) if (kappa .ne. kappaBoundNz) kappaHarm = kappaHarm*BR - B(I) = B(I) + (kappaHarm) * T_Bathz2 + B(I) = B(I) + (kappaHarm) * T_Bathz2 !+ boundray_term_vel(ix,iy,nz,T_Bathz2) end if end if end if @@ -180,5 +180,49 @@ function constantboundarytempgrad(I) result(temp) ! so we need to add T_BathCG to the temperature at the boundary end function constantboundarytempgrad !!!############################################################################################### + + !!!######################################################################## + !!! This subroutine calculates the value of the convective term of the H matrix... + !!! ...at the boundary. + !!!######################################################################## + function boundray_term_vel(x_b, y_b, z_b, TB) result(vel_conv) + integer(int12), intent(in) :: x_b, y_b, z_b + real(real12) :: vel_conv + real(real12) :: rho, CV, TB + real(real12), dimension(3) :: vel_in, vel_out + + !------------------------------------------------------------ + ! The boundary term is calculated of the boundary grid point. + !------------------------------------------------------------ + + + rho = grid(x_b,y_b,z_b)%rho + CV = grid(x_b,y_b,z_b)%heat_capacity + vel_in = grid(x_b,y_b,z_b)%vel + vel_out = grid(x_b,y_b,z_b)%vel + + if (x_b .eq. 1_int12) then + vel_conv = vel_in(1)*rho*CV*(1.0_real12/(2.0_real12*grid(x_b,y_b,z_b)%Length(1))) + vel_conv = -vel_conv + else if (x_b .eq. nx) then + vel_conv = vel_in(1)*rho*CV*(1.0_real12/(2.0_real12*grid(x_b,y_b,z_b)%Length(1))) + + else if (y_b .eq. 1_int12) then + vel_conv = vel_in(2)*rho*CV*(1.0_real12/(2.0_real12*grid(x_b,y_b,z_b)%Length(2))) + vel_conv = -vel_conv + else if (y_b .eq. ny) then + vel_conv = vel_in(2)*rho*CV*(1.0_real12/(2.0_real12*grid(x_b,y_b,z_b)%Length(2))) + else if (z_b .eq. 1_int12) then + vel_conv = vel_in(3)*rho*CV*(1.0_real12/(2.0_real12*grid(x_b,y_b,z_b)%Length(3))) + vel_conv = -vel_conv + else if (z_b .eq. nz) then + vel_conv = vel_in(3)*rho*CV*(1.0_real12/(2.0_real12*grid(x_b,y_b,z_b)%Length(3))) + else + vel_conv = 0.0_real12 + end if + vel_conv = vel_conv*TB + !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + end function boundray_term_vel + !!!######################################################################## end module boundary_vector diff --git a/src/heatflow/mod_constructions.f90 b/src/heatflow/mod_constructions.f90 index b233f7c..6e8afc0 100644 --- a/src/heatflow/mod_constructions.f90 +++ b/src/heatflow/mod_constructions.f90 @@ -26,19 +26,21 @@ module constructions integer(int12) :: imaterial_type !what type of material it is in the mat.in real(real12) :: volume !volume of the block real(real12), dimension(3) :: Length !length of the block in 3 dimensions + real(real12), dimension(3) :: vel !vel of the block in 3 dimensions integer(int12) :: iheater !whether the block is a heater - real(real12) :: kappa, rho, heat_capacity, tau, em !physical properties of the material + real(real12) :: kappa, rho, heat_capacity, tau, em !physical properties of the material end type heatblock !!Defines the material to be used in the simulation type material integer(int12) :: index !identifier for the material real(real12) :: heat_capacity !heat capacity of the material - real(real12) :: h_conv !convective heat transfer coefficient + ! real(real12) :: h_conv !convective heat transfer coefficient real(real12) :: kappa !thermal conductivity - real(real12) :: kappa3D !three-dimensional thermal conductivity + ! real(real12) :: kappa3D !three-dimensional thermal conductivity real(real12) :: rho !density of the material - real(real12) :: sound_speed !speed of sound in the material + ! real(real12) :: sound_speed !speed of sound in the material real(real12) :: tau !relaxation time + real(real12), dimension(3) :: vel !velocity of the material real(real12) :: em !emisstivity logical :: source !whether the material is a source of heat ?? end type material diff --git a/src/heatflow/mod_evolve.f90 b/src/heatflow/mod_evolve.f90 index 1209b25..0554c0d 100644 --- a/src/heatflow/mod_evolve.f90 +++ b/src/heatflow/mod_evolve.f90 @@ -26,14 +26,20 @@ module evolution use sptype, only: I4B use solver, only: linbcg use globe_data, only: Temp_p, Temp_pp, inverse_time, heat, lin_rhoc + use globe_data, only: acsr, ja, ia use heating, only: heater use boundary_vector, only: boundary use cattaneo, only: S_catS - use tempdep, only: ChangeProp +! use tempdep, only: ChangeProp + use petsc_solver, only: solve_petsc_csr + implicit none private public :: simulate + + ! Module-level variables for PETSc (persist across time steps) + integer, allocatable, save :: ia32(:), ja32(:) ! 32-bit copies for PETSc contains @@ -44,11 +50,13 @@ module evolution !!!################################################################################################# subroutine simulate(itime) integer(int12), intent(in) :: itime - real(real12), dimension(NA) :: S, x, Q, Qdens, S_CAT, B - integer(int12) :: ncg, itol, itmax !, iss - integer(I4B) :: iter + real(real12), dimension(NA) :: S, Q, Qdens, S_CAT, B + real(real12), dimension(:), allocatable :: x + integer:: ncg, itol, itmax !, iss + integer :: iter real(real12) :: e, err, tol - + integer :: NA32 + !---------------------- ! Initialize vectors !---------------------- @@ -111,6 +119,16 @@ subroutine simulate(itime) !--------------------------------------------- if ( iSteady .eq. 0 ) then S = - inverse_time * Temp_p * lin_rhoc - Qdens - B + if (IVERB .gt. 3) then + write(*,*) "S construction diagnostics:" + write(*,*) " inverse_time =", inverse_time + write(*,*) " Temp_p avg =", sum(Temp_p)/size(Temp_p) + write(*,*) " lin_rhoc avg =", sum(lin_rhoc)/size(lin_rhoc) + write(*,*) " Qdens avg =", sum(Qdens)/size(Qdens) + write(*,*) " B avg =", sum(B)/size(B) + write(*,*) " -inverse_time*Temp_p*lin_rhoc avg =", sum(-inverse_time*Temp_p*lin_rhoc)/size(Temp_p) + write(*,*) " S before S_CAT avg =", sum(S)/size(S) + end if if ( iCAttaneo .eq. 1) then S = S + S_CAT end if @@ -137,19 +155,59 @@ subroutine simulate(itime) ! iter: Output - gives the number of the final iteration. ! err: Output - records the error of the final iteration. ! iss: Input - sets the Sparse Storage type (1=SRS, 2=SDS). - x=Temp_p+(Temp_p-Temp_pp) - if (any(x-Temp_p .lt. TINY)) x=x+TINY !avoid nan solver issue + ! x=Temp_p+(Temp_p-Temp_pp) + ! if (any(x-Temp_p .lt. TINY)) x=x+TINY !avoid nan solver issue itol=1 tol=1.e-32_real12 - itmax=50000 + itmax=500000 ncg = 0 - iter=ncg + iter= 0 err=E - - CALL linbcg(S,x,itol=int(itol,I4B),tol=tol, itmax=int(itmax,I4B), iter=iter, & - err=E) + ! call bicgstab(acsr, ia, ja, S, itmax, Temp_p, x, iter) + + ! Allocate and initialize x with a good initial guess + allocate(x(NA)) + x = Temp_p + (Temp_p - Temp_pp) + if (any(x - Temp_p .lt. TINY)) x = x + TINY ! avoid nan solver issue + + ! Debug: Print initial guess statistics + if (IVERB .gt. 3) then + write(*,*) "========== PETSc Solver Diagnostics ==========" + write(*,*) "Time step:", itime + write(*,*) "Initial guess x: min=", minval(x), " max=", maxval(x), " avg=", sum(x)/size(x) + write(*,*) "RHS S: min=", minval(S), " max=", maxval(S), " avg=", sum(S)/size(S) + write(*,*) "Temp_p: min=", minval(Temp_p), " max=", maxval(Temp_p), " avg=", sum(Temp_p)/size(Temp_p) + write(*,*) "Matrix acsr: min=", minval(acsr), " max=", maxval(acsr), " avg=", sum(acsr)/size(acsr) + write(*,*) "Matrix size: n=", NA32, " nnz=", size(acsr) + end if + + ! Convert to 32-bit integers for PETSc (only on first call) + if (.not. allocated(ia32)) then + allocate(ia32(size(ia)), ja32(size(ja))) + ia32 = int(ia, kind=kind(ia32)) + ja32 = int(ja, kind=kind(ja32)) + end if + NA32 = int(NA, kind=kind(NA32)) + + call solve_petsc_csr(NA32, ia32, ja32, acsr, S, x, tol, itmax) + + ! Debug: Print solution statistics and verify solution + if (IVERB .gt. 3) then + write(*,*) "Solution x after PETSc: min=", minval(x), " max=", maxval(x), " avg=", sum(x)/size(x) + write(*,*) "Temperature change: avg(x-Temp_p)=", sum(x-Temp_p)/size(x) + write(*,*) "Max temperature change: ", maxval(abs(x-Temp_p)) + write(*,*) "==============================================" + end if + + ! Note: Don't deallocate ia32, ja32 - keep them for next time step + + + ! CALL solve_pardiso(acsr, S, ia, ja, x) + ! CALL linbcg(S,x,itol=int(itol,I4B),tol=tol, itmax=int(itmax,I4B), iter=iter, & + ! err=E) + ! if (any(isnan(x(:)))) then write(0,*) "fatal error: NAN in x tempurature vector" write(0,*) 'time step ', itime, " T ", sum(Temp_p)/size(Temp_p), E ,iter @@ -169,9 +227,9 @@ subroutine simulate(itime) Temp_pp = Temp_p Temp_p = x - if (TempDepProp .eq. 1) then - CALL ChangeProp() - end if + ! if (TempDepProp .eq. 1) then + ! CALL ChangeProp() + ! end if !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ end subroutine simulate diff --git a/src/heatflow/mod_global.f90 b/src/heatflow/mod_global.f90 index 5207578..439952c 100644 --- a/src/heatflow/mod_global.f90 +++ b/src/heatflow/mod_global.f90 @@ -4,7 +4,7 @@ !!!################################################################################################# module globe_data use sptype, only: sprs2_dp, diag_sprs_dp - use constants, only: real12 + use constants, only: real12, int12 implicit none real(real12), dimension(:,:,:), allocatable :: Temp_cur ! Current temperature real(real12), dimension(:), allocatable :: Temp_p, Temp_pp ! Previous and previous previous temperature @@ -13,6 +13,10 @@ module globe_data real(real12), dimension(:), allocatable :: lin_rhoc ! 1D array for HeatCapacity*Rho TYPE(sprs2_dp) :: ra !Techniqually rH TYPE(diag_sprs_dp) :: da !Techniqually dH + real(real12), dimension(:), allocatable :: acsr + integer(kind=8), dimension(:), allocatable :: ja + integer(kind=8), dimension(:), allocatable :: ia + character(len=1024) :: logname ! Log file name end module globe_data diff --git a/src/heatflow/mod_heating.f90 b/src/heatflow/mod_heating.f90 index 831396f..8fb0298 100644 --- a/src/heatflow/mod_heating.f90 +++ b/src/heatflow/mod_heating.f90 @@ -14,7 +14,7 @@ module Heating use constants, only: real12, int12, pi, StefBoltz use globe_data, only: Temp_p, Temp_pp, Heat, heated_volume, Q_P, heated_temp use inputs, only: nx,ny,nz, grid, NA, power_in, time_step, heated_steps, T_System, freq, ntime, & - T_Bath + T_Bath, icattaneo use materials, only: material implicit none contains @@ -30,7 +30,7 @@ subroutine heater(itime, Q, Qdens) integer(int12) :: ix, iy, iz, IA ,heated_num real(real12) :: time, POWER, time_pulse, x, x2 real(real12) :: rho, volume, heat_capacity, area, tau, sum_temp - + logical :: shared_power ! Initialize variables IA = 0 Q = 0._real12 @@ -53,7 +53,7 @@ subroutine heater(itime, Q, Qdens) heat_capacity = grid(ix,iy,iz)%heat_capacity area = grid(ix,iy,iz)%Length(1)*grid(ix,iy,iz)%Length(2) !??? !tau divided by time_step squared in setup.f90 - tau = grid(ix,iy,iz)%tau*(time_step**2.0_real12) + tau = grid(ix,iy,iz)%tau*(time_step) ! select heater case select case(grid(ix,iy,iz)%iheater) @@ -78,9 +78,24 @@ subroutine heater(itime, Q, Qdens) !------------------------------ if ( time .le. time_pulse ) then Q(IA) = POWER + ! print *, "Heating on" else + ! print *, "Heating off" Q(IA) = 0.0_real12 end if + + if (icattaneo .eq. 1) then + if (itime .eq. 1) then + Q(IA) = Q(IA) + tau*POWER + ! print *, "Turning on heater at time ", time + end if + if (itime .eq. (heated_steps + 1)) then + Q(IA) = Q(IA) - tau*POWER + ! print *, "Turning off heater at time ", time + end if + end if + + !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ case(3) !------------------------------ @@ -114,13 +129,50 @@ subroutine heater(itime, Q, Qdens) Q(IA) = 0.0 end if !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + case(7) + !--------------------------------------------------------- + ! Heater on for a time period, off for a time period, + ! then on again (square-wave heating) + !--------------------------------------------------------- + + if (mod(time, 2.0_real12 * time_pulse) .le. time_pulse) then + Q(IA) = POWER + + if (icattaneo .eq. 1) then + + ! Forward-time pulse check + if (mod(time + time_step, 2.0_real12 * time_pulse) .le. time_pulse) then + Q(IA) = Q(IA) + tau * POWER + end if + + ! Backward-time pulse check + if (mod(time - time_step, 2.0_real12 * time_pulse) .le. time_pulse) then + Q(IA) = Q(IA) - tau * POWER + end if + + end if + + else + Q(IA) = 0.0_real12 + end if + + + case(10) + Q(IA) = POWER - (tau*(POWER)) + + case(11) + Q(IA) = POWER + + case(12) + Q(IA) = POWER + (tau*(POWER)) end select !------------------------------ ! If emissitivity is not zero, then calculate the radiative heating !------------------------------ - Q(IA) = Q(IA) - grid(ix,iy,iz)%em * grid(ix,iy,iz)%length(1)*& - grid(ix,iy,iz)%length(2)*StefBoltz & - * ((Temp_p(IA)**4.0_real12) - (T_Bath**4.0_real12)) + ! Q(IA) = Q(IA) - grid(ix,iy,iz)%em * grid(ix,iy,iz)%length(1)*& + ! grid(ix,iy,iz)%length(2)*StefBoltz & + ! * ((Temp_p(IA)**4.0_real12) - (T_Bath**4.0_real12)) !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -161,11 +213,15 @@ subroutine heater(itime, Q, Qdens) ! Normalize all heat sources by the heated volume - if (heated_volume .gt. 0.0) then - Qdens(:) = Q(:) / heated_volume - heated_temp = sum_temp / heated_volume - end if - + shared_power = .TRUE. + if (shared_power) then + if (heated_volume .gt. 0.0) then + Qdens(:) = Q(:) / heated_volume + heated_temp = sum_temp / heated_volume + end if + else + Qdens(:) = Q(:) / grid(1,1,1)%volume + end if end subroutine heater diff --git a/src/heatflow/mod_hmatrix.f90 b/src/heatflow/mod_hmatrix.f90 index 9e0b219..fb7b747 100644 --- a/src/heatflow/mod_hmatrix.f90 +++ b/src/heatflow/mod_hmatrix.f90 @@ -148,7 +148,8 @@ function hmatrixfunc(i, j) result(H) if (x .eq. 1) then H=0.0_real12 else - H = A ! X left neighbor (left cell interaction) + H = A ! X left neighbor (left cell interaction) + !H = H + calculate_convective_conductivity(xm, y, z, x, y, z) end if end if @@ -157,6 +158,7 @@ function hmatrixfunc(i, j) result(H) H=0.0_real12 else H = B ! X right neighbor (right cell interaction) + !H = H + calculate_convective_conductivity(xp, y, z, x, y, z) end if end if @@ -165,6 +167,7 @@ function hmatrixfunc(i, j) result(H) H=0.0_real12 else H = D ! Y down neighbor (down cell interaction) + !H = H + calculate_convective_conductivity(x, ym, z, x, y, z) end if end if if ((i-j) .eq. -nx) then @@ -172,6 +175,7 @@ function hmatrixfunc(i, j) result(H) H=0.0_real12 else H = E ! Y up neighbor (up cell interaction) + !H = H + calculate_convective_conductivity(x, yp, z, x, y, z) end if end if @@ -181,6 +185,7 @@ function hmatrixfunc(i, j) result(H) else !write(*,*) 'F this is forward (in) z',F H = F ! Z in neighbor (forward cell interaction) !!!Frank had this as G during testing + !H = H + calculate_convective_conductivity(x, y, zm, x, y, z) end if end if @@ -190,6 +195,7 @@ function hmatrixfunc(i, j) result(H) else !write(*,*) 'G this is backward (out) z?',G H = G ! Z out neighbor (backward cell interaction) !!!Frank Had this as F during testing + !H = H + calculate_convective_conductivity(x, y, zp, x, y, z) end if end if @@ -391,5 +397,118 @@ subroutine boundry_diag_term(x_b, y_b, z_b, x, y, z, kappa_ab) end subroutine boundry_diag_term !!!######################################################################## + !!!######################################################################## + !!! subroutine to calculate the convective conductivity between two points + !!!######################################################################## + function calculate_convective_conductivity(x_in, y_in, z_in, x_out, y_out, z_out) result(vel_conv) + integer(int12), intent(in) :: x_in, y_in, z_in, x_out, y_out, z_out + real(real12) :: vel_conv + real(real12) :: rho, CV + real(real12), dimension(3) :: vel_in, vel_out + + vel_conv = 0.0_real12 + + ! if not an edge element + if ((x_in .ge. 1) .and. (x_in .le. nx) .and. (y_in .ge. 1) .and. & + (y_in .le. ny) .and. (z_in .ge. 1) .and. (z_in .le. nz)) then + + rho = grid(x_out,y_out,z_out)%rho + CV = grid(x_out,y_out,z_out)%heat_capacity + vel_in = grid(x_in,y_in,z_in)%vel + vel_out = grid(x_out,y_out,z_out)%vel + + if (x_in .ne. x_out) then + if (vel_in(1) .ne. vel_out(1)) then + vel_conv = 0.0_real12 + else + vel_conv = vel_in(1)*rho*CV*(1.0_real12/(2.0_real12*grid(x_out,y_out,z_out)%Length(1))) + if (x_in .lt. x_out) then + vel_conv = -vel_conv + end if + end if + + else if (y_in .ne. y_out) then + if (vel_in(2) .ne. vel_out(2)) then + vel_conv = 0.0_real12 + else + vel_conv = vel_in(2)*rho*CV*(1.0_real12/(2.0_real12*grid(x_out,y_out,z_out)%Length(2))) + if (y_in .lt. y_out) then + vel_conv = -vel_conv + end if + end if + + else if (z_in .ne. z_out) then + if (vel_in(3) .ne. vel_out(3)) then + vel_conv = 0.0_real12 + else + vel_conv = vel_in(3)*rho*CV*(1.0_real12/(2.0_real12*grid(x_out,y_out,z_out)%Length(3))) + if (z_in .lt. z_out) then + vel_conv = -vel_conv + end if + end if + end if + + else + vel_conv = 0.0_real12 + end if + vel_conv = vel_conv + end function calculate_convective_conductivity + + !!!######################################################################## + + !!!######################################################################## + !!! This subroutine calculates the value of the convective term of the H matrix... + !!! ...at the boundary. + !!!######################################################################## + subroutine boundry_diag_term_vel(x_b, y_b, z_b, x, y, z, vel_conv) + integer(int12), intent(in) :: x_b, y_b, z_b, x, y, z + real(real12), intent(out) :: vel_conv + real(real12) :: rho, CV + real(real12), dimension(3) :: vel_in, vel_out + + !------------------------------------------------------------ + ! The boundary term is calculated of the boundary grid point. + !------------------------------------------------------------ + + + rho = grid(x,y,z)%rho + CV = grid(x,y,z)%heat_capacity + vel_in = grid(x_b,y_b,z_b)%vel + vel_out = grid(x,y,z)%vel + + if (x_b .eq. 1) then + if (vel_in(1) .ne. vel_out(1)) then + vel_conv = 0.0_real12 + else + vel_conv = vel_in(1)*rho*CV*(1.0_real12/(2.0_real12*grid(x,y,z)%Length(1))) + if (x_b .lt. x) then + vel_conv = -vel_conv + end if + end if + + else if (y_b .ne. y) then + if (vel_in(2) .ne. vel_out(2)) then + vel_conv = 0.0_real12 + else + vel_conv = vel_in(2)*rho*CV*(1.0_real12/(2.0_real12*grid(x,y,z)%Length(2))) + if (y_b .lt. y) then + vel_conv = -vel_conv + end if + end if + + else if (z_b .ne. z) then + if (vel_in(3) .ne. vel_out(3)) then + vel_conv = 0.0_real12 + else + vel_conv = vel_in(3)*rho*CV*(1.0_real12/(2.0_real12*grid(x,y,z)%Length(3))) + if (z_b .lt. z) then + vel_conv = -vel_conv + end if + end if + end if + !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + end subroutine boundry_diag_term_vel + !!!######################################################################## + end module hmatrixmod diff --git a/src/heatflow/mod_inputs.f90 b/src/heatflow/mod_inputs.f90 index 4d29245..7f1fd9b 100644 --- a/src/heatflow/mod_inputs.f90 +++ b/src/heatflow/mod_inputs.f90 @@ -75,6 +75,7 @@ module inputs use constants, only: real12, int12 use constructions, only: heatblock, material + use mpi implicit none integer :: unit, newunit @@ -91,7 +92,7 @@ module inputs integer(int12) :: start_ix, end_ix, start_iy, end_iy, start_iz, end_iz, TempDepProp, heated_steps ! flags logical :: Check_Sparse_Full, Check_Stability, Check_Steady_State - logical :: WriteToTxt, LPercentage, InputTempDis + logical :: WriteToTxt, LPercentage, InputTempDis, CompressedOutput logical :: Test_Run = .FALSE., FullRestart = .FALSE. ! Name of simiulation run @@ -204,7 +205,7 @@ end subroutine read_all_files subroutine read_param(unit) implicit none integer:: unit, Reason - integer,dimension(45)::readvar + integer,dimension(46)::readvar character(1024)::buffer readvar(:)=0 @@ -221,6 +222,7 @@ subroutine read_param(unit) RunName = 'default' RunName = trim(adjustl(RunName)) WriteToTxt = .FALSE. + CompressedOutput = .FALSE. ntime = 10 heated_steps = 0 write_every = 1 @@ -325,6 +327,7 @@ subroutine read_param(unit) CALL assignD(buffer,"BR",BR,readvar(43)) CALL assignL(buffer,"T_BathCC",T_BathCC,readvar(44)) CALL assignS(buffer,"CG_dir",CG_dir,readvar(45)) + CALL assignL(buffer,"_CompressedOutput",CompressedOutput,readvar(46)) !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ end do @@ -365,18 +368,24 @@ end subroutine read_param !!!################################################################################################# subroutine check_param(readvar,n) implicit none - integer::n,i + integer::n,i,ierr,comm_rank integer,dimension(n)::readvar + logical :: is_root ! Not currently in use + call MPI_Comm_rank(MPI_COMM_WORLD, comm_rank, ierr) + is_root = (comm_rank == 0) + if(any(readvar.gt.1)) then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## ERROR ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- Error in subroutine "check_param" ---' - write(6,'(A)') ' --- ERROR: same KEYWORD apears more than once ---' + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## ERROR ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- Error in subroutine "check_param" ---' + write(6,'(A)') ' --- ERROR: same KEYWORD apears more than once ---' + end if stop end if @@ -387,13 +396,15 @@ subroutine check_param(readvar,n) !------------------------------------------------------------------------------------ ErrKB:if (((any(readvar(9:11).eq.0)) .or. any(readvar(29:31).eq.0)) & .and. (readvar(28) .eq. 0) )then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## ERORR ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- Error in subroutine "check_param" ---' - write(6,'(A)') ' --- ERROR: KappaBoundx,y,z and KappaBound are not set ---' + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## ERORR ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- Error in subroutine "check_param" ---' + write(6,'(A)') ' --- ERROR: KappaBoundx,y,z and KappaBound are not set ---' + end if stop end if ErrKB !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -403,13 +414,15 @@ subroutine check_param(readvar,n) WarKBO:if (((all(readvar(9:11).eq.1)) .or. all(readvar(29:31).eq.1)) & .and. (readvar(28) .eq. 0) )then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## Warning ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- Warning in subroutine "check_param" ---' - write(6,'(A)') ' --- Warning: KappaBound is not set ---' + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## Warning ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- Warning in subroutine "check_param" ---' + write(6,'(A)') ' --- Warning: KappaBound is not set ---' + end if readvar(28) = 1 end if WarKBO @@ -419,13 +432,15 @@ subroutine check_param(readvar,n) !------------------------------------------------------------------------------------ WarKB:if (((any(readvar(9:11).eq.0)) .or. any(readvar(29:31).eq.0)) & .and. (readvar(28) .eq. 1) )then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## WARNING ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- Warning in subroutine "check_param" ---' - write(6,'(A)') ' --- WARNING: KappaBoundx,y,z not set, using KappaBound ---' + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## WARNING ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- Warning in subroutine "check_param" ---' + write(6,'(A)') ' --- WARNING: KappaBoundx,y,z not set, using KappaBound ---' + end if kappaBoundx1 = KappaBound kappaBoundy1 = KappaBound kappaBoundz1 = KappaBound @@ -436,39 +451,45 @@ subroutine check_param(readvar,n) !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ elseif (((any(readvar(9:11).eq.0)) .or. any(readvar(29:31).eq.0)) & .or. (readvar(28) .gt. 0) ) then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## WARNING ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- Warning in subroutine "check_param" ---' - write(6,'(A)') ' --- WARNING: Periodic Boundry set, set Kappa are ignored ---' + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## WARNING ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- Warning in subroutine "check_param" ---' + write(6,'(A)') ' --- WARNING: Periodic Boundry set, set Kappa are ignored ---' + end if end if PB !------------------------------------------------------------------------------------ ! warning about missing bath temps. reassine to T_Bath !------------------------------------------------------------------------------------ if ((readvar(42) .eq. 1) .and. (T_BathCG .gt. 0)) then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## WARNING ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- Warning in subroutine "check_param" ---' - write(6,'(A)') ' --- WARNING: T_BathCG set T_Bath/ T_Bath x,y,z will not be used ---' + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## WARNING ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- Warning in subroutine "check_param" ---' + write(6,'(A)') ' --- WARNING: T_BathCG set T_Bath/ T_Bath x,y,z will not be used ---' + end if !set all T_Bath value checks to 1 readvar(20:25) = 1 readvar(27) = 1 end if WarBath:if ( any(readvar(20:25).eq.0) ) then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## WARNING ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- Warning in subroutine "check_param" ---' - write(6,'(A)') ' --- WARNING: T_Bath x,y,z not set T_Bath will be used ---' + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## WARNING ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- Warning in subroutine "check_param" ---' + write(6,'(A)') ' --- WARNING: T_Bath x,y,z not set T_Bath will be used ---' + end if T_Bathx1 = T_Bath T_Bathx2 = T_Bath T_Bathy1 = T_Bath @@ -479,13 +500,15 @@ subroutine check_param(readvar,n) !Check if T_BathCG less than 0 if ((readvar(42) .eq. 1) .and. (T_BathCG .lt. 0.0_real12)) then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## WARNING ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- Warning in subroutine "check_param" ---' - write(6,'(A)') ' --- WARNING: T_BathCG is negative, are you sure you? ---' + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## WARNING ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- Warning in subroutine "check_param" ---' + write(6,'(A)') ' --- WARNING: T_BathCG is negative, are you sure you? ---' + end if end if if (readvar(42) .eq. 0) then T_BathCG = 0 @@ -498,43 +521,52 @@ subroutine check_param(readvar,n) ! Further warnings !------------------------------------------------------------------------------------ if (any(readvar(32:37).eq.0)) then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## WARNING ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- WARNING in subroutine "check_param" ---' - write(6,'(A)') ' --- WARNING: Some or All output write cells paramters are not defined ---' - write(6,*) ' --- USING: ', 'Start_ix = ', start_ix, ', end_ix = ', end_ix, ', start_iy = ', & + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## WARNING ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- WARNING in subroutine "check_param" ---' + write(6,'(A)') ' --- WARNING: Some or All output write cells paramters are not defined ---' + write(6,*) ' --- USING: ', 'Start_ix = ', start_ix, ', end_ix = ', end_ix, ', start_iy = ', & start_iy,', end_iy = ', end_iy, ', start_iz = ', start_iz, ', end_iz = ', end_iz + end if end if if (readvar(39) .eq. 0) then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## WARNING ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- WARNING in subroutine "check_param" ---' - write(6,'(A)') ' --- WARNING: TempDepProp not set, no action needed ---' + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## WARNING ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- WARNING in subroutine "check_param" ---' + write(6,'(A)') ' --- WARNING: TempDepProp not set, no action needed ---' + end if readvar(39) = 1 end if + if (readvar(46) .eq. 0) then + readvar(46) = 1 + end if + if (any(readvar.eq.0)) then - write(6,*) - write(6,'(A43)') '###############################' - write(6,'(A43)') '########## WARNING ##########' - write(6,'(A43)') '###############################' - write(6,*) - write(6,'(A)') ' --- WARNING in subroutine "check_param" ---' - write(6,'(A)') ' --- WARNING: Essential parameters missing ---' - ! Print all indices of readvar that are 0 - do i = 1, size(readvar) - if (readvar(i) == 0) then + if (is_root) then + write(6,*) + write(6,'(A43)') '###############################' + write(6,'(A43)') '########## WARNING ##########' + write(6,'(A43)') '###############################' + write(6,*) + write(6,'(A)') ' --- WARNING in subroutine "check_param" ---' + write(6,'(A)') ' --- WARNING: Essential parameters missing ---' + do i = 1, size(readvar) + if (readvar(i) == 0) then write(6, '(A,I3)') 'Index ', i, ' of readvar is 0' - end if - end do - write(6,*) + end if + end do + write(6,*) + end if end if !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -542,7 +574,7 @@ subroutine check_param(readvar,n) !------------------------------------------------------------------------------------ ! verbos to check for errors !------------------------------------------------------------------------------------ - if(IVERB .ge. 1) then + if(is_root .and. IVERB .ge. 1) then write(6,'(A)') ' vebose printing option' write(6,'(A)') ' running calculation with :' write(6,'(A35,I6)') ' IVERB = ', IVERB @@ -554,6 +586,7 @@ subroutine check_param(readvar,n) write(6,'(A35,L1)') ' _FullRestart = ', FullRestart write(6,'(A35,A)') ' _RunName = ', trim(RunName) write(6,'(A35,L1)') ' _WriteToTxt = ', WriteToTxt + write(6,'(A35,L1)') ' _CompressedOutput = ', CompressedOutput write(6,'(A35,I12)') ' ntime = ', ntime write(6,'(A35,I12)') ' heattime = ', heated_steps write(6,'(A35,I12)') ' write_every = ', write_every @@ -657,7 +690,7 @@ subroutine read_mat(unit) type(material), dimension(100) :: dum_mat character(1024) :: buffer integer :: reason, j - integer, dimension(8) :: readvarmat + integer, dimension(6) :: readvarmat integer :: i, index i=0 @@ -697,13 +730,11 @@ subroutine read_mat(unit) CALL assignD(buffer,"heat_capacity",dum_mat(i)%heat_capacity,readvarmat(1))! assign heatCapacity - CALL assignD(buffer,"h_conv" ,dum_mat(i)%h_conv ,readvarmat(2))! assign h_conv - CALL assignD(buffer,"kappa" ,dum_mat(i)%kappa ,readvarmat(3))! assign kappa - CALL assignD(buffer,"kappa3D" ,dum_mat(i)%kappa3D ,readvarmat(4))! assign kappa3D - CALL assignD(buffer,"rho" ,dum_mat(i)%rho ,readvarmat(5))! assign rho - CALL assignD(buffer,"sound_speed" ,dum_mat(i)%sound_speed ,readvarmat(6))! assign sound_speed - CALL assignD(buffer,"tau" ,dum_mat(i)%tau ,readvarmat(7))! assign tau - CALL assignD(buffer,"em" ,dum_mat(i)%em ,readvarmat(8))! assign e + CALL assignD(buffer,"kappa" ,dum_mat(i)%kappa ,readvarmat(2))! assign kappa + CALL assignD(buffer,"rho" ,dum_mat(i)%rho ,readvarmat(3))! assign rho + CALL assignD(buffer,"tau" ,dum_mat(i)%tau ,readvarmat(4))! assign tau + CALL assignD(buffer,"em" ,dum_mat(i)%em ,readvarmat(5))! assign e + CALL assignV(buffer,"vel" ,dum_mat(i)%vel ,readvarmat(6)) ! assign velocity end do read ! Check for duplicate indices @@ -827,5 +858,23 @@ function val(buffer) end function val !!!################################################################################################# +!!!################################################################################################# +!!! assign velocity +!!!################################################################################################# + subroutine assignV(buffer, keyword, variable, found) + implicit none + integer::found + character(1024)::buffer1,buffer2 + character(*)::buffer,keyword + real(real12), dimension(3)::variable + buffer1=buffer(:scan(buffer,"=")-1) + if(scan("=",buffer).ne.0) buffer2=val(buffer) + if(trim(adjustl(buffer1)).eq.trim(adjustl(keyword))& + .and.trim(adjustl(buffer2)).ne.'') then + found=found+1 + read(buffer2,*) variable(1), variable(2), variable(3) + end if + end subroutine assignV + end module inputs diff --git a/src/heatflow/mod_material.f90 b/src/heatflow/mod_material.f90 index 6ca7192..65b9fdf 100644 --- a/src/heatflow/mod_material.f90 +++ b/src/heatflow/mod_material.f90 @@ -51,12 +51,13 @@ module materials !!! - em, the emissivity of the material. !!!######################################################################### -subroutine material(imaterial_type,kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau, em) +subroutine material(imaterial_type,kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau, em, vel) integer(int12), intent(in) ::imaterial_type integer(int12) :: i, tmp real(real12), intent(inout) :: kappa3D, kappa, h_conv, heat_capacity, sound_speed, rho, tau real(real12), intent(inout) :: em + real(real12), dimension(3), intent(inout) :: vel logical :: found !!!============================================= @@ -90,12 +91,15 @@ subroutine material(imaterial_type,kappa,kappa3D,h_conv,heat_capacity,rho,sound_ if (tmp .eq. imaterial_type) then found = .true. heat_capacity = input_materials(i)%heat_capacity - h_conv = input_materials(i)%h_conv + ! h_conv = input_materials(i)%h_conv kappa = input_materials(i)%kappa - kappa3D = input_materials(i)%kappa3D + ! kappa3D = input_materials(i)%kappa3D rho = input_materials(i)%rho - sound_speed = input_materials(i)%sound_speed + ! sound_speed = input_materials(i)%sound_speed tau = input_materials(i)%tau + vel(1) = input_materials(i)%vel(1) + vel(2) = input_materials(i)%vel(2) + vel(3) = input_materials(i)%vel(3) em = input_materials(i)%em exit mat_loop diff --git a/src/heatflow/mod_output.f90 b/src/heatflow/mod_output.f90 index edafa9a..c518761 100644 --- a/src/heatflow/mod_output.f90 +++ b/src/heatflow/mod_output.f90 @@ -44,9 +44,9 @@ module output use constants, only: real12, int12, TINY, fields use inputs, only: nx,ny,nz, time_step, grid, NA, Check_Steady_State, ntime, WriteToTxt - use inputs, only: Test_Run, freq, RunName, FullRestart, IVERB, write_every + use inputs, only: Test_Run, freq, RunName, FullRestart, IVERB, write_every, CompressedOutput use inputs, only: start_ix, end_ix, start_iy, end_iy, start_iz, end_iz - use globe_data, only: Temp_p,Temp_pp, heat, heated_volume + use globe_data, only: Temp_p,Temp_pp, heat, heated_volume, logname implicit none contains @@ -54,9 +54,11 @@ subroutine data_write(itime) implicit none integer(int12), intent(in) :: itime real(real12), dimension(nx,ny,nz) :: CT, Temp_cur - integer(int12) :: ix, iy, iz, indexA, logunit - character(len=1024) :: file_prefix, file_extension, outdir, logname + integer(int12) :: ix, iy, iz, indexA + character(len=1024) :: file_prefix, file_extension, outdir + integer :: logunit + logunit = 20 file_prefix = 'Temperture_' outdir='./outputs/' file_extension = '.out' @@ -90,9 +92,18 @@ subroutine data_write(itime) ! find most recent log file and open it !--------------------------------------- CALL last_log(logname,outdir) - open(newunit=logunit,file=logname) - write(logunit,*) real((itime-1)*(time_step)), & - (Temp_cur(start_ix:end_ix, start_iy:end_iy, start_iz:end_iz)) + if (CompressedOutput) then + open(logunit,file=logname, status='unknown', access='stream', position='append') + write(logunit) real((itime-1)*(time_step), kind=real12) + write(logunit) (Temp_cur(start_ix:end_ix, start_iy:end_iy, start_iz:end_iz)) + else + open(logunit,file=logname) + ! print*, logunit + ! print*, logname + write(logunit,*) real((itime-1)*(time_step)), & + (Temp_cur(start_ix:end_ix, start_iy:end_iy, start_iz:end_iz)) + end if + close(logunit) !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ end if end if @@ -105,13 +116,23 @@ subroutine data_write(itime) ! write out to log file !--------------------------------------- if (.not. Test_run) then - if (WriteToTxt) then - if (mod(itime, write_every) .eq. 0) then - write(*, *) 'Writing Temperature difference to file' - write(logunit,*) real((itime-1)*(time_step)), & - (Temp_cur(start_ix:end_ix, start_iy:end_iy, start_iz:end_iz)) - end if - endif + if (mod(itime, write_every) .eq. 0) then + if (CompressedOutput) then + write(*, *) 'Writing Temperature (Compressed) to file' + open(logunit,file=logname, status='unknown', access='stream', position='append') + write(logunit) real((itime-1)*(time_step), kind=real12) + write(logunit) (Temp_cur(start_ix:end_ix, start_iy:end_iy, start_iz:end_iz)) + close(logunit) + elseif (WriteToTxt) then + write(*, *) 'Writing Temperature difference to file' + ! print*, logunit + ! print*, logname + open(logunit,file=logname, status='old', position='append') + write(logunit,*) real((itime-1)*(time_step)), & + (Temp_cur(start_ix:end_ix, start_iy:end_iy, start_iz:end_iz)) + close(logunit) + end if + end if end if !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -193,8 +214,13 @@ subroutine last_log(logname,outdir) i = 0 flag=.true. do while (flag) - write(logname, '(A,A,I2.2)') trim(adjustl(outdir)) // 'output_' // & - trim(adjustl(RunName)),'_', i + if (CompressedOutput) then + write(logname, '(A,A,I2.2,A)') trim(adjustl(outdir)) // 'output_' // & + trim(adjustl(RunName)),'_', i, '.bin' + else + write(logname, '(A,A,I2.2)') trim(adjustl(outdir)) // 'output_' // & + trim(adjustl(RunName)),'_', i + endif inquire(file=logname, exist=flag) i = i+1 end do diff --git a/src/heatflow/mod_petsc_solver.f90 b/src/heatflow/mod_petsc_solver.f90 new file mode 100644 index 0000000..65bbf45 --- /dev/null +++ b/src/heatflow/mod_petsc_solver.f90 @@ -0,0 +1,337 @@ +module petsc_solver +#include "petsc/finclude/petscsys.h" +#include "petsc/finclude/petscksp.h" + use petscksp + use mpi + implicit none + private + public :: petsc_init, petsc_finalize, solve_petsc_csr, petsc_cleanup, petsc_is_root, petsc_world_size + + ! ===== PRECONDITIONER SELECTION ===== + ! Change this to switch between preconditioners: + ! 'GAMG' = Algebraic Multigrid (best for elliptic PDEs, 10-20x faster) + ! 'ILU' = Incomplete LU (good general purpose, robust) + ! 'LU' = Direct solver (most robust, uses more memory) + character(len=10), parameter :: PRECONDITIONER = 'GAMG' + ! ==================================== + + ! Persistent PETSc objects (reused across timesteps for memory efficiency) + Mat, save :: A_saved = PETSC_NULL_MAT + Vec, save :: bb_saved = PETSC_NULL_VEC + Vec, save :: xx_saved = PETSC_NULL_VEC + KSP, save :: ksp_saved = PETSC_NULL_KSP + logical, save :: initialized = .false. + integer, save :: n_saved = 0 + integer, save :: comm_rank_saved = 0 + integer, save :: comm_size_saved = 1 + integer, save :: row_start_saved = 1 + integer, save :: row_end_saved = 0 + PetscInt, allocatable, save :: ia_saved(:), ja_saved(:) + PetscInt, allocatable, save :: local_indices_saved(:) + PetscScalar, allocatable, save :: aval_saved(:) + PetscInt, allocatable, save :: diag_nnz_saved(:), offdiag_nnz_saved(:) + integer, allocatable, save :: recvcounts_saved(:), displs_saved(:) + real(8), allocatable, save :: b_local_saved(:), x_local_saved(:) + +contains + + subroutine petsc_init() + integer :: ierr + call PetscInitialize(PETSC_NULL_CHARACTER, ierr) + call MPI_Comm_rank(PETSC_COMM_WORLD, comm_rank_saved, ierr) + call MPI_Comm_size(PETSC_COMM_WORLD, comm_size_saved, ierr) + end subroutine petsc_init + + logical function petsc_is_root() + petsc_is_root = (comm_rank_saved == 0) + end function petsc_is_root + + integer function petsc_world_size() + petsc_world_size = comm_size_saved + end function petsc_world_size + + subroutine petsc_finalize() + integer :: ierr + call petsc_cleanup() + call PetscFinalize(ierr) + end subroutine petsc_finalize + + subroutine petsc_cleanup() + integer :: ierr + + if (A_saved /= PETSC_NULL_MAT) call MatDestroy(A_saved, ierr) + if (bb_saved /= PETSC_NULL_VEC) call VecDestroy(bb_saved, ierr) + if (xx_saved /= PETSC_NULL_VEC) call VecDestroy(xx_saved, ierr) + if (ksp_saved /= PETSC_NULL_KSP) call KSPDestroy(ksp_saved, ierr) + if (allocated(ia_saved)) deallocate(ia_saved) + if (allocated(ja_saved)) deallocate(ja_saved) + if (allocated(local_indices_saved)) deallocate(local_indices_saved) + if (allocated(aval_saved)) deallocate(aval_saved) + if (allocated(diag_nnz_saved)) deallocate(diag_nnz_saved) + if (allocated(offdiag_nnz_saved)) deallocate(offdiag_nnz_saved) + if (allocated(recvcounts_saved)) deallocate(recvcounts_saved) + if (allocated(displs_saved)) deallocate(displs_saved) + if (allocated(b_local_saved)) deallocate(b_local_saved) + if (allocated(x_local_saved)) deallocate(x_local_saved) + + A_saved = PETSC_NULL_MAT + bb_saved = PETSC_NULL_VEC + xx_saved = PETSC_NULL_VEC + ksp_saved = PETSC_NULL_KSP + initialized = .false. + n_saved = 0 + row_start_saved = 1 + row_end_saved = 0 + end subroutine petsc_cleanup + + subroutine compute_partition(n, rank, nproc, row_start, row_end) + integer, intent(in) :: n, rank, nproc + integer, intent(out) :: row_start, row_end + integer :: base_rows, remainder_rows, local_rows + + base_rows = n / nproc + remainder_rows = mod(n, nproc) + local_rows = base_rows + if (rank < remainder_rows) local_rows = local_rows + 1 + + row_start = rank * base_rows + min(rank, remainder_rows) + 1 + row_end = row_start + local_rows - 1 + end subroutine compute_partition + + subroutine preallocate_local_rows(n, ia, ja) + integer, intent(in) :: n + integer, intent(in) :: ia(:), ja(:) + integer :: local_row, global_row, entry_idx, nlocal + integer :: diag_begin, diag_end + + nlocal = max(0, row_end_saved - row_start_saved + 1) + allocate(diag_nnz_saved(nlocal), offdiag_nnz_saved(nlocal)) + diag_nnz_saved = 0 + offdiag_nnz_saved = 0 + + diag_begin = row_start_saved + diag_end = row_end_saved + do local_row = 1, nlocal + global_row = row_start_saved + local_row - 1 + do entry_idx = ia(global_row), ia(global_row + 1) - 1 + if (ja(entry_idx) >= diag_begin .and. ja(entry_idx) <= diag_end) then + diag_nnz_saved(local_row) = diag_nnz_saved(local_row) + 1 + else + offdiag_nnz_saved(local_row) = offdiag_nnz_saved(local_row) + 1 + end if + end do + end do + end subroutine preallocate_local_rows + + subroutine build_distributed_matrix(n, ia, ja, aval) + integer, intent(in) :: n + integer, intent(in) :: ia(:), ja(:) + real(8), intent(in) :: aval(:) + + PetscInt :: row_idx(1) + PetscInt, allocatable :: cols0(:) + PetscScalar, allocatable :: vals0(:) + integer :: ierr, global_row, local_row, nlocal, row_nnz, max_row_nnz + + call preallocate_local_rows(n, ia, ja) + nlocal = max(0, row_end_saved - row_start_saved + 1) + max_row_nnz = max(1, maxval(ia(2:n + 1) - ia(1:n))) + + call MatCreate(PETSC_COMM_WORLD, A_saved, ierr) + call MatSetSizes(A_saved, nlocal, nlocal, n, n, ierr) + call MatSetType(A_saved, MATAIJ, ierr) + call MatSeqAIJSetPreallocation(A_saved, 0, diag_nnz_saved, ierr) + call MatMPIAIJSetPreallocation(A_saved, 0, diag_nnz_saved, 0, offdiag_nnz_saved, ierr) + + allocate(cols0(max_row_nnz), vals0(max_row_nnz)) + do local_row = 1, nlocal + global_row = row_start_saved + local_row - 1 + row_nnz = ia(global_row + 1) - ia(global_row) + if (row_nnz <= 0) cycle + + row_idx(1) = global_row - 1 + cols0(1:row_nnz) = ja(ia(global_row):ia(global_row + 1) - 1) - 1 + vals0(1:row_nnz) = aval(ia(global_row):ia(global_row + 1) - 1) + call MatSetValues(A_saved, 1, row_idx, row_nnz, cols0, vals0, INSERT_VALUES, ierr) + end do + deallocate(cols0, vals0) + + call MatAssemblyBegin(A_saved, MAT_FINAL_ASSEMBLY, ierr) + call MatAssemblyEnd(A_saved, MAT_FINAL_ASSEMBLY, ierr) + end subroutine build_distributed_matrix + + subroutine update_distributed_matrix(n, ia, ja, aval) + integer, intent(in) :: n + integer, intent(in) :: ia(:), ja(:) + real(8), intent(in) :: aval(:) + + PetscInt :: row_idx(1) + PetscInt, allocatable :: cols0(:) + PetscScalar, allocatable :: vals0(:) + integer :: ierr, global_row, local_row, nlocal, row_nnz, max_row_nnz + + nlocal = max(0, row_end_saved - row_start_saved + 1) + max_row_nnz = max(1, maxval(ia(2:n + 1) - ia(1:n))) + allocate(cols0(max_row_nnz), vals0(max_row_nnz)) + + call MatZeroEntries(A_saved, ierr) + do local_row = 1, nlocal + global_row = row_start_saved + local_row - 1 + row_nnz = ia(global_row + 1) - ia(global_row) + if (row_nnz <= 0) cycle + + row_idx(1) = global_row - 1 + cols0(1:row_nnz) = ja(ia(global_row):ia(global_row + 1) - 1) - 1 + vals0(1:row_nnz) = aval(ia(global_row):ia(global_row + 1) - 1) + call MatSetValues(A_saved, 1, row_idx, row_nnz, cols0, vals0, INSERT_VALUES, ierr) + end do + deallocate(cols0, vals0) + + call MatAssemblyBegin(A_saved, MAT_FINAL_ASSEMBLY, ierr) + call MatAssemblyEnd(A_saved, MAT_FINAL_ASSEMBLY, ierr) + end subroutine update_distributed_matrix + + subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit) + integer, intent(in) :: n + integer, intent(in) :: ia(:), ja(:) + real(8), intent(in) :: aval(:), b(:) + real(8), intent(inout) :: x(:) + real(8), intent(in) :: rtol + integer, intent(in) :: maxit + + PC :: pc + integer :: ierr, its + real(8) :: rnorm + logical :: rebuild_needed + integer :: nlocal, rank_idx, local_row + + if (size(ia) /= n + 1) stop 'solve_petsc_csr: ia size mismatch' + if (size(b) /= n .or. size(x) /= n) stop 'solve_petsc_csr: vector size mismatch' + + rebuild_needed = .false. + if (.not. initialized) rebuild_needed = .true. + if (n /= n_saved) rebuild_needed = .true. + + if (.not. rebuild_needed .and. allocated(ia_saved)) then + if (size(ia_saved) /= size(ia) .or. size(ja_saved) /= size(ja) .or. size(aval_saved) /= size(aval)) then + rebuild_needed = .true. + else if (.not. all(ia_saved == ia - 1) .or. .not. all(ja_saved == ja - 1)) then + rebuild_needed = .true. + end if + end if + + if (rebuild_needed) then + if (initialized) call petsc_cleanup() + + call compute_partition(n, comm_rank_saved, comm_size_saved, row_start_saved, row_end_saved) + nlocal = max(0, row_end_saved - row_start_saved + 1) + + allocate(ia_saved(size(ia)), ja_saved(size(ja)), aval_saved(size(aval))) + ia_saved = ia - 1 + ja_saved = ja - 1 + aval_saved = aval + + allocate(local_indices_saved(max(1, nlocal))) + if (nlocal > 0) then + do local_row = 1, nlocal + local_indices_saved(local_row) = row_start_saved + local_row - 2 + end do + end if + + call build_distributed_matrix(n, ia, ja, aval) + + allocate(recvcounts_saved(comm_size_saved), displs_saved(comm_size_saved)) + do rank_idx = 0, comm_size_saved - 1 + call compute_partition(n, rank_idx, comm_size_saved, ierr, its) + recvcounts_saved(rank_idx + 1) = max(0, its - ierr + 1) + displs_saved(rank_idx + 1) = ierr - 1 + end do + + allocate(b_local_saved(max(1, nlocal)), x_local_saved(max(1, nlocal))) + + call VecCreateMPI(PETSC_COMM_WORLD, nlocal, n, bb_saved, ierr) + call VecDuplicate(bb_saved, xx_saved, ierr) + + call KSPCreate(PETSC_COMM_WORLD, ksp_saved, ierr) + call KSPSetOperators(ksp_saved, A_saved, A_saved, ierr) + call KSPSetInitialGuessNonzero(ksp_saved, PETSC_TRUE, ierr) + call KSPGetPC(ksp_saved, pc, ierr) + + select case (trim(PRECONDITIONER)) + case ('GAMG') + call PCSetType(pc, PCGAMG, ierr) + call KSPSetType(ksp_saved, KSPGMRES, ierr) + if (petsc_is_root()) then + write(*,'(A,I0,A)') ' [Solver] Using GAMG (Algebraic Multigrid) with GMRES across ', & + comm_size_saved, ' MPI ranks' + end if + + case ('ILU') + call PCSetType(pc, PCILU, ierr) + call KSPSetType(ksp_saved, KSPBCGS, ierr) + if (petsc_is_root()) then + write(*,'(A,I0,A)') ' [Solver] Using ILU preconditioner with BiCGSTAB across ', & + comm_size_saved, ' MPI ranks' + end if + + case ('LU') + call PCSetType(pc, PCLU, ierr) + call KSPSetType(ksp_saved, KSPPREONLY, ierr) + if (petsc_is_root()) then + write(*,'(A,I0,A)') ' [Solver] Using direct LU solver across ', comm_size_saved, ' MPI ranks' + end if + + case default + if (petsc_is_root()) then + write(*,'(A,A)') ' [Warning] Unknown preconditioner: ', trim(PRECONDITIONER) + write(*,'(A)') ' Defaulting to ILU' + end if + call PCSetType(pc, PCILU, ierr) + call KSPSetType(ksp_saved, KSPBCGS, ierr) + end select + + call KSPSetTolerances(ksp_saved, rtol, PETSC_DEFAULT_REAL, PETSC_DEFAULT_REAL, maxit, ierr) + call KSPSetNormType(ksp_saved, KSP_NORM_UNPRECONDITIONED, ierr) + call KSPSetFromOptions(ksp_saved, ierr) + + initialized = .true. + n_saved = n + else + if (any(aval_saved /= aval)) then + aval_saved = aval + call update_distributed_matrix(n, ia, ja, aval) + end if + end if + + nlocal = max(0, row_end_saved - row_start_saved + 1) + if (nlocal > 0) then + b_local_saved(1:nlocal) = b(row_start_saved:row_end_saved) + x_local_saved(1:nlocal) = x(row_start_saved:row_end_saved) + end if + + call VecSet(bb_saved, 0.0d0, ierr) + call VecSet(xx_saved, 0.0d0, ierr) + if (nlocal > 0) then + call VecSetValues(bb_saved, nlocal, local_indices_saved, b_local_saved, INSERT_VALUES, ierr) + call VecSetValues(xx_saved, nlocal, local_indices_saved, x_local_saved, INSERT_VALUES, ierr) + end if + call VecAssemblyBegin(bb_saved, ierr) + call VecAssemblyEnd(bb_saved, ierr) + call VecAssemblyBegin(xx_saved, ierr) + call VecAssemblyEnd(xx_saved, ierr) + + call KSPSolve(ksp_saved, bb_saved, xx_saved, ierr) + if (ierr /= 0) then + write(0,*) 'ERROR: KSPSolve failed with error code:', ierr + stop + end if + + call KSPGetIterationNumber(ksp_saved, its, ierr) + call KSPGetResidualNorm(ksp_saved, rnorm, ierr) + + if (nlocal > 0) call VecGetValues(xx_saved, nlocal, local_indices_saved, x_local_saved, ierr) + + call MPI_Allgatherv(x_local_saved, nlocal, MPI_DOUBLE_PRECISION, x, recvcounts_saved, displs_saved, & + MPI_DOUBLE_PRECISION, PETSC_COMM_WORLD, ierr) + end subroutine solve_petsc_csr +end module petsc_solver \ No newline at end of file diff --git a/src/heatflow/mod_setup.f90 b/src/heatflow/mod_setup.f90 index 17496a9..57cd2f7 100644 --- a/src/heatflow/mod_setup.f90 +++ b/src/heatflow/mod_setup.f90 @@ -16,15 +16,19 @@ !!!################################################################################################# module setup use constants, only: real12, int12, TINY + use mpi use inputs, only: nx, ny, nz, NA, grid, time_step, kappaBoundx1, kappaBoundy1, kappaBoundz1 use inputs, only: Check_Sparse_Full, Check_Stability, ntime,IVERB, Periodicx, Periodicy use inputs, only: Periodicz ! use hmatrixmod, only: hmatrixfunc use globe_data, only: ra, Temp_cur, Temp_p, Temp_pp,inverse_time, heat, lin_rhoc, Q_P + use globe_data, only: acsr, ja, ia use solver, only: SRSin use materials, only: material implicit none + public :: set_global_variables + contains @@ -33,7 +37,13 @@ module setup !!!################################################################################################# subroutine set_global_variables() integer(int12) :: ix,iy,iz,index + integer :: ierr, comm_rank + logical :: is_root real(real12) :: kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau, em + real(real12), dimension(3) :: vel + + call MPI_Comm_rank(MPI_COMM_WORLD, comm_rank, ierr) + is_root = (comm_rank == 0) allocate(Temp_cur(nx, ny, nz)) allocate(Temp_p(NA)) @@ -48,19 +58,27 @@ subroutine set_global_variables() ! can be expanded to include more properties at a ! later date !--------------------------------------------------- - write(*,*) "Setting up material properties" + if (is_root) then + write(*,*) "Setting up material properties" + write(*,'(A,I10,A)') " Processing ", NA, " grid cells..." + end if index = 0 do iz = 1, nz + ! Progress reporting every 10% for large grids + if (is_root .and. mod(iz-1, max(1,nz/10)) == 0 .and. iz > 1) then + write(*,'(A,I3,A)') " Progress: ", int(100.0*real(iz)/real(nz)), "%" + end if do iy = 1, ny do ix = 1, nx index = index + 1 CALL material(grid(ix,iy,iz)%imaterial_type,& - kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau, em) + kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau, em, vel) grid(ix,iy,iz)%kappa = kappa grid(ix,iy,iz)%rho = rho grid(ix,iy,iz)%heat_capacity = heat_capacity grid(ix,iy,iz)%tau = tau*inverse_time*inverse_time grid(ix,iy,iz)%em = em + grid(ix,iy,iz)%vel(:) = vel(:) lin_rhoc(index) = rho*heat_capacity if (Check_Stability) CALL stability(kappa, rho, heat_capacity, ix, iy, iz) end do @@ -71,11 +89,15 @@ subroutine set_global_variables() !--------------------------------------------------- ! Check if the sparse matrix matches the full matrix !--------------------------------------------------- + if (is_root) write(*,*) "Building sparse H matrix..." if (Check_Sparse_Full) then print*, "CHECK SPARSE FULL" CALL build_Hmatrix() else + ! Build CSR format directly (acsr, ja, ia are allocated inside sparse_Hmatrix) CALL sparse_Hmatrix() + ! No need for COO->CSR conversion anymore, it's already in CSR format! + if (is_root) write(*,*) "Sparse matrix setup complete." end if !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -84,26 +106,28 @@ end subroutine set_global_variables !!!################################################################################################# !!!################################################################################################# -!!! This sets up the H Matrix directly in sparse row storage +!!! This sets up the H Matrix directly in sparse row storage (CSR format) +!!! Modified to build CSR directly instead of COO->CSR to save memory !!!################################################################################################# subroutine sparse_Hmatrix() implicit none + integer :: ierr, comm_rank + logical :: is_root real(real12) :: H0 ! Holds the value of the H matrix - integer(int12) :: i, j, len, count, k ! i and j are the row and column of the H matrix + integer(int12) :: i, j, count, k, row ! i and j are the row and column of the H matrix + integer(int12) :: nnz_estimate ! Holds the values to add to the row to get the column integer(int12), allocatable, dimension(:) :: addit - ! The number of non-zero elements in the H matrix to look for - len = 7*nx*ny*nz - 2*(nx*ny + ny*nz + nz*nx) - if (Periodicx) len = len + 2*ny*nz - if (Periodicy) len = len + 2*nz*nx - if (Periodicz) len = len + 2*nx*ny + ! Temporary arrays for building each row + real(real12), allocatable, dimension(:) :: row_vals + integer(int12), allocatable, dimension(:) :: row_cols + integer(int12) :: row_count, max_row_size + + call MPI_Comm_rank(MPI_COMM_WORLD, comm_rank, ierr) + is_root = (comm_rank == 0) + ra%n = NA ! The number of rows in the H matrix - ra%len = len ! The number of non-zero elements in the H matrix - ! Allocate the arrays to hold the H matrix in sparse storage - allocate(ra%val(len), ra%irow(len), ra%jcol(len)) - ra%val(:)=0 - ra%irow(:)=-2 - ra%jcol(:)=-1 + ra%len = 0 addit = [1] ! The values to add to the row to get the column if (Periodicx) addit = [addit, (nx-1)] if (ny .gt. 1) addit = [addit, nx] ! Add the values to add to the row to get the column @@ -111,47 +135,149 @@ subroutine sparse_Hmatrix() if (nz .gt. 1) addit = [addit, nx*ny] ! Add the values to add to the row to get the column if ((Periodicz).and.(nz .gt. 1)) addit = [addit, (nz-1)*ny*nx] + max_row_size = 1 + 2*size(addit,1) + nnz_estimate = NA * max_row_size + + if (allocated(acsr)) deallocate(acsr) + if (allocated(ja)) deallocate(ja) + if (allocated(ia)) deallocate(ia) + allocate(acsr(nnz_estimate), ja(nnz_estimate), ia(NA+1)) + allocate(row_vals(max_row_size), row_cols(max_row_size)) + acsr = 0.0_real12 + ja = 0 + ia = 0 + !write(6,*) NA, nx,ny,nz !write(6,*) addit,size(addit,1) !write(6,*) NA !write(6,*) "=========================================" - count = 0 ! The number of non-zero elements in the H matrix - parent_loop: do j = 1, NA ! Loop over the columns of the H matrix - i=j ! The row of the H matrix - count = count + 1 ! The number of non-zero elements in the H matrix - H0 = hmatrixfunc(i,j) ! The value of the H matrix - ra%val(count) = H0 ! The value of the H matrix - ra%irow(count) = i ! The row of the H matrix - ra%jcol(count) = j ! The column of the H matrix - ! Loop over the values to add to the row to get the column - neighbour_loop: do k = 1, size(addit,1) - i = j + addit(k) ! The row of the H matrix - ! If the row is greater than the number of rows ... - !...in the H matrix then go to the next column - if ((i.gt.NA)) cycle parent_loop - H0=hmatrixfunc(i,j) ! The value of the H matrix - ! If the value of the H matrix is less than TINY then go to the next value ... - !...to add to the row to get the column - if (abs(H0).lt.TINY) cycle neighbour_loop - count = count + 1 ! The number of non-zero elements in the H matrix - ra%val(count) = H0 ! The value of the H matrix - ra%irow(count) = i ! The row of the H matrix - ra%jcol(count) = j ! The column of the H matrix - count = count + 1 ! The number of non-zero elements in the H matrix - ra%val(count) = H0 ! The value of the H matrix - ra%irow(count) = j ! The row of the H matrix - ra%jcol(count) = i ! The column of the H matrix - !write(6,*) j,i, H0, count - end do neighbour_loop - end do parent_loop - - !write(6,*) "=========================================" - !write(6,*) 'c len',count, len + + count = 0 ! Total nonzeros counter + ia(1) = 1 ! CSR row pointer (1-based for Fortran) + + ! Build CSR format row-by-row + if (is_root) write(*,'(A)') " Building CSR matrix row-by-row..." + parent_loop: do row = 1, NA + ! Progress reporting every 10% + if (is_root .and. mod(row-1, max(1,NA/10)) == 0 .and. row > 1) then + write(*,'(A,I3,A,I12,A)') " Progress: ", int(100.0*real(row)/real(NA)), & + "%, nnz=", count, "" + end if + + row_count = 0 + + ! Diagonal element + j = row + row_count = row_count + 1 + H0 = hmatrixfunc(row, j) + row_vals(row_count) = H0 + row_cols(row_count) = j + + ! Off-diagonal elements (process in column-sorted order for CSR) + ! First pass: collect all neighbors + do k = 1, size(addit,1) + j = row + addit(k) + if (j > NA) cycle ! Skip if out of bounds + + H0 = hmatrixfunc(row, j) + if (abs(H0) >= TINY) then + row_count = row_count + 1 + row_vals(row_count) = H0 + row_cols(row_count) = j + end if + end do + + ! Second pass: collect reverse neighbors (j < row) + do k = 1, size(addit,1) + j = row - addit(k) + if (j < 1) cycle ! Skip if out of bounds + + H0 = hmatrixfunc(row, j) + if (abs(H0) >= TINY) then + row_count = row_count + 1 + row_vals(row_count) = H0 + row_cols(row_count) = j + end if + end do + + ! Sort this row's entries by column index (required for CSR) + call sort_row(row_vals, row_cols, row_count) + + ! Copy row data to CSR arrays (no bounds checking - we pre-allocated correctly) + do i = 1, row_count + count = count + 1 + acsr(count) = row_vals(i) + ja(count) = row_cols(i) + end do + + ! Update row pointer + ia(row+1) = count + 1 + end do parent_loop + + ra%len = count + + ! Trim arrays to actual size if we over-estimated + if (count < size(acsr)) then + if (is_root) write(*,'(A,I12,A,I12)') " Trimming arrays from ", size(acsr), " to ", count + call trim_csr_arrays(acsr, ja, count) + end if + + deallocate(addit, row_vals, row_cols) + if (is_root) write(*,'(A,I12,A)') " CSR matrix built successfully. Actual nonzeros: ", count, "" end subroutine sparse_Hmatrix !!!################################################################################################# +!!!################################################################################################# +!!! Sort a row's entries by column index (simple insertion sort, rows are small) +!!!################################################################################################# + subroutine sort_row(vals, cols, n) + implicit none + integer(int12), intent(in) :: n + real(real12), dimension(n), intent(inout) :: vals + integer(int12), dimension(n), intent(inout) :: cols + integer(int12) :: i, j, temp_col + real(real12) :: temp_val + + do i = 2, n + temp_val = vals(i) + temp_col = cols(i) + j = i - 1 + do while (j >= 1) + if (cols(j) <= temp_col) exit + vals(j+1) = vals(j) + cols(j+1) = cols(j) + j = j - 1 + end do + vals(j+1) = temp_val + cols(j+1) = temp_col + end do + end subroutine sort_row +!!!################################################################################################# + +!!!################################################################################################# +!!! Trim CSR arrays to exact size +!!!################################################################################################# + subroutine trim_csr_arrays(acsr_arr, ja_arr, final_size) + implicit none + integer(int12), intent(in) :: final_size + real(real12), allocatable, dimension(:), intent(inout) :: acsr_arr + integer(int12), allocatable, dimension(:), intent(inout) :: ja_arr + real(real12), allocatable, dimension(:) :: temp_vals + integer(int12), allocatable, dimension(:) :: temp_cols + + allocate(temp_vals(final_size), temp_cols(final_size)) + temp_vals = acsr_arr(1:final_size) + temp_cols = ja_arr(1:final_size) + deallocate(acsr_arr, ja_arr) + allocate(acsr_arr(final_size), ja_arr(final_size)) + acsr_arr = temp_vals + ja_arr = temp_cols + deallocate(temp_vals, temp_cols) + end subroutine trim_csr_arrays +!!!################################################################################################# + !!!################################################################################################# !!! Check if the simulation will be stable. NOT FULLY IMPLEMENTED !!!################################################################################################# diff --git a/src/heatflow/mod_tempdep.f90 b/src/heatflow/mod_tempdep.f90 old mode 100755 new mode 100644 index 475617c..b8b936a --- a/src/heatflow/mod_tempdep.f90 +++ b/src/heatflow/mod_tempdep.f90 @@ -27,7 +27,7 @@ module TempDep use inputs, only: Grid, TempDepProp, Nz, Ny, Nx - use setup, only: sparse_Hmatrix + ! use setup, only: sparse_Hmatrix use globe_data, only: Temp_p, lin_rhoc use constants, only: real12, int12 @@ -35,34 +35,34 @@ module TempDep contains - subroutine ChangeProp() - character(len=100) :: filename - integer(int12) :: ix,iy,iz, index - logical :: res - index = 1 - !Loop over all the grid points - do iz = 1, Nz - do iy = 1, Ny - do ix = 1, Nx - ! Construct the filename for the Material table asscoiated with the grid point - filename = trim('./inputs/MatTable' // & - trim(adjustl(char(Grid(ix, iy, iz)%imaterial_type)))) - inquire(file=filename, exist=res) - if (res) then - ! Read the temperature dependent properties from the file - ! CALL ReadTempDepTable(filename, ix, iy, iz, index) - else - ! File does not exist, continue to the next grid point - continue - end if - ! Read the temperature dependent properties from the file - index = index + 1 - end do - end do - end do - ! Construct the sparse matrix - CALL sparse_Hmatrix() - end subroutine ChangeProp + ! subroutine ChangeProp() + ! character(len=100) :: filename + ! integer(int12) :: ix,iy,iz, index + ! logical :: res + ! index = 1 + ! !Loop over all the grid points + ! do iz = 1, Nz + ! do iy = 1, Ny + ! do ix = 1, Nx + ! ! Construct the filename for the Material table asscoiated with the grid point + ! filename = trim('./inputs/MatTable' // & + ! trim(adjustl(char(Grid(ix, iy, iz)%imaterial_type)))) + ! inquire(file=filename, exist=res) + ! if (res) then + ! ! Read the temperature dependent properties from the file + ! ! CALL ReadTempDepTable(filename, ix, iy, iz, index) + ! else + ! ! File does not exist, continue to the next grid point + ! continue + ! end if + ! ! Read the temperature dependent properties from the file + ! index = index + 1 + ! end do + ! end do + ! end do + ! ! Construct the sparse matrix + ! CALL sparse_Hmatrix() + ! end subroutine ChangeProp ! subroutine ReadTempDepTable(filename, ix, iy, iz, index) ! character(len=*), intent(in) :: filename