diff --git a/.gitignore b/.gitignore
index b4431b4..259e90f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,6 +6,8 @@ cases/
 .#*
 *.mod
 build/
+build-*/
+build_*/
 inputs/
 outputs/
 DElectrodes/
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 31395de..70f5058 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,56 +1,93 @@
 cmake_minimum_required(VERSION 3.17.5...3.27.7)
 
-# define build environments
+# ========================
+# Project metadata
+# ========================
 set(PROJECT_NAME "HeatFlow")
 set(PROJECT_LIBRARY_NAME "HeatFlow_library")
 
-set( CMAKE_INSTALL_PREFIX "$ENV{HOME}/.local/${PROJECT_NAME}"
-  CACHE STRING "Select where to install the program." )
+project(HeatFlow
+  DESCRIPTION "Fortran heat flow"
+  LANGUAGES C Fortran)
+
+set(CMAKE_INSTALL_PREFIX "$ENV{HOME}/.local/${PROJECT_NAME}"
+  CACHE STRING "Install prefix")
+
 set(CMAKE_BUILD_DIR ${CMAKE_CURRENT_BINARY_DIR}
-  CACHE STRING "Select where to build the program." )
-set(MODULE_DIR ${CMAKE_BUILD_DIR}/mod)
+  CACHE STRING "Build directory")
 
+set(MODULE_DIR ${CMAKE_BUILD_DIR}/mod)
 
-# set compiler
-set(CMAKE_Fortran_COMPILER gfortran
-  CACHE STRING "Select Fortran compiler." )  # Change this to your desired compiler
-set(CMAKE_C_COMPILER gcc
-  CACHE STRING "Select C compiler." )  # Change this to your desired compiler
 set(CMAKE_Fortran_STANDARD 2018)
 
-# set the project name
-project(HeatFlow LANGUAGES C Fortran)
-
-# set options for building tests and examples
-option(BUILD_TESTS "Build the unit tests" On)
+option(BUILD_TESTS "Build unit tests" ON)
+
+# ========================
+# Build types
+# ========================
+set(CMAKE_CONFIGURATION_TYPES
+  "Release" "Parallel" "Serial" "Dev" "Debug" "Coverage" "Bigmem"
+  CACHE STRING "Build types")
+
+set(CMAKE_BUILD_TYPE "Release"
+  CACHE STRING "Default build type")
+
+# ========================
+# Compiler detection
+# ========================
+if(NOT DEFINED CMAKE_Fortran_COMPILER)
+  find_program(DEFAULT_Fortran_COMPILER NAMES mpifort mpif90 gfortran ifx ifort nagfor)
+  if(DEFAULT_Fortran_COMPILER)
+    set(CMAKE_Fortran_COMPILER "${DEFAULT_Fortran_COMPILER}"
+      CACHE STRING "Fortran compiler")
+  endif()
+endif()
 
-# set coverage compiler flags
-if (CMAKE_BUILD_TYPE MATCHES "Coverage")
-  list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake")
-  set(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake)
-  if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
-      include(CodeCoverage)
-      setup_target_for_coverage_gcovr_html(
-        NAME coverage
-        EXECUTABLE ctest
-        EXCLUDE "${CMAKE_SOURCE_DIR}/test/*")
+if(NOT DEFINED CMAKE_C_COMPILER)
+  find_program(DEFAULT_C_COMPILER NAMES mpicc gcc clang icx icc)
+  if(DEFAULT_C_COMPILER)
+    set(CMAKE_C_COMPILER "${DEFAULT_C_COMPILER}"
+      CACHE STRING "C compiler")
   endif()
 endif()
 
-# enable testing
-enable_testing()
+# ========================
+# Flags
+# ========================
+if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
+    set(PPFLAGS -cpp)
+    set(MPFLAGS -fopenmp)
+    set(WARNFLAGS -Wall)
+    set(DEVFLAGS -g -fbacktrace -fcheck=all -fbounds-check -Og)
+    set(DEBUGFLAGS -fbounds-check)
+    set(MEMFLAGS -mcmodel=large)
+    set(OPTIMFLAGS -O3 -march=native)
 
-set( PROJECT_DESCRIPTION
-  "Fortran heat flow" )
-set( PROJECT_URL "https://github.com/ExeQuantCode/HeatFlow" )
-set( CMAKE_CONFIGURATION_TYPES "Release" "Parallel" "Serial" "Dev" "Debug" "Coverage" "Bigmem"
-  CACHE STRING "List of configurations types." )
-set( CMAKE_BUILD_TYPE "Release"
-  CACHE STRING "Select which configuration to build." )
+elseif (CMAKE_Fortran_COMPILER_ID STREQUAL "NAG")
+    set(PPFLAGS -f2018 -fpp)
+    set(MPFLAGS -openmp)
+    set(WARNFLAGS -Wall)
+    set(DEVFLAGS -g -C=all -O0)
+    set(DEBUGFLAGS -C=array)
+    set(MEMFLAGS -mcmodel=large)
+    set(OPTIMFLAGS -O3)
+
+elseif (CMAKE_Fortran_COMPILER_ID MATCHES "Intel")
+    set(PPFLAGS -fpp)
+    set(MPFLAGS -qopenmp)
+    set(WARNFLAGS -warn all)
+    set(DEVFLAGS -check all -warn)
+    set(DEBUGFLAGS -check all -fpe0 -traceback)
+    set(MEMFLAGS -mcmodel=large)
+    set(OPTIMFLAGS -O3)
+endif()
 
-# Define the sources
+# ========================
+# Sources
+# ========================
 set(SRC_DIR src)
 set(LIB_DIR ${SRC_DIR}/heatflow)
+set(PETSC_DIR_SRC ${SRC_DIR}/petsc_solver)
 
 set(LIB_FILES
     mod_constants.f90
@@ -76,157 +113,174 @@ set(SRC_FILES
 )
 
 foreach(lib ${LIB_FILES})
-    list(APPEND PREPENDED_LIB_FILES ${CMAKE_CURRENT_LIST_DIR}/${LIB_DIR}/${lib})
+  list(APPEND PREPENDED_LIB_FILES ${CMAKE_CURRENT_LIST_DIR}/${LIB_DIR}/${lib})
 endforeach()
 
 foreach(src ${SRC_FILES})
-    list(APPEND PREPENDED_SRC_FILES ${CMAKE_CURRENT_LIST_DIR}/${SRC_DIR}/${src})
+  list(APPEND PREPENDED_SRC_FILES ${CMAKE_CURRENT_LIST_DIR}/${SRC_DIR}/${src})
 endforeach()
 
+# PETSc solver now isolated
+set(PETSC_SOLVER_FILE
+  ${CMAKE_CURRENT_LIST_DIR}/${PETSC_DIR_SRC}/mod_petsc_solver.f90)
+set(PETSC_SOLVER_OBJECT ${CMAKE_CURRENT_BINARY_DIR}/mod_petsc_solver.o)
 
-# initialise flags
-set(CPPFLAGS "")
-set(CFLAGS "")
-set(MODULEFLAGS "")
-set(MPFLAGS "")
-set(WARNFLAGS "")
-set(DEVFLAGS "")
-set(DEBUGFLAGS "")
-set(MEMFLAGS "")
-set(OPTIMFLAGS "")
-set(FASTFLAGS "")
-
-# set flags based on compiler
-if (CMAKE_Fortran_COMPILER MATCHES ".*gfortran.*" OR CMAKE_Fortran_COMPILER MATCHES ".*gcc.*")
-    message(STATUS "Using gfortran compiler")
-    set(PPFLAGS -cpp)
-    set(MPFLAGS -fopenmp)
-    set(WARNFLAGS -Wall)
-    set(DEVFLAGS -g -fbacktrace -fcheck=all -fbounds-check -Og)
-    set(DEBUGFLAGS -fbounds-check)
-    set(MEMFLAGS -mcmodel=large)
-    set(OPTIMFLAGS -O3 -march=native)
-    set(FASTFLAGS -Ofast -march=native)
-elseif (CMAKE_Fortran_COMPILER MATCHES ".*nag.*")
-    message(STATUS "Using nag compiler")
-    set(PPFLAGS -f2018 -fpp)
-    set(MPFLAGS -openmp)
-    set(WARNFLAGS -Wall)
-    set(DEVFLAGS -g -mtrace -C=all -colour -O0)
-    set(DEBUGFLAGS -C=array)
-    set(MEMFLAGS -mcmodel=large)
-    set(OPTIMFLAGS -O3)
-    set(FASTFLAGS -Ofast)
-elseif (CMAKE_Fortran_COMPILER MATCHES ".*ifort.*" OR CMAKE_Fortran_COMPILER MATCHES ".*ifx.*")
-    message(STATUS "Using intel compiler")
-    set(PPFLAGS -fpp)
-    set(MPFLAGS -qopenmp)
-    set(WARNFLAGS -warn all)
-    set(DEVFLAGS -check all -warn)
-    set(DEBUGFLAGS -check all -fpe0 -warn -tracekback -debug extended)
-    set(MEMFLAGS -mcmodel=large)
-    set(OPTIMFLAGS -O3)
-    set(FASTFLAGS -Ofast)
+# ========================
+# MPI
+# ========================
+get_filename_component(FORTRAN_COMPILER_NAME "${CMAKE_Fortran_COMPILER}" NAME)
+
+add_library(HeatFlowMPI INTERFACE)
+add_library(HeatFlow::MPI ALIAS HeatFlowMPI)
+
+if(FORTRAN_COMPILER_NAME MATCHES "^mpi(fort|f90)$")
+  message(STATUS "Using MPI wrapper ${CMAKE_Fortran_COMPILER}")
 else()
-    # Code for other Fortran compilers
-    message(STATUS "Using a different Fortran compiler")
+  find_package(MPI REQUIRED COMPONENTS Fortran)
+  target_link_libraries(HeatFlowMPI INTERFACE MPI::MPI_Fortran)
 endif()
 
+set(PETSC_MPI_INCLUDE_FLAGS "")
+if(NOT FORTRAN_COMPILER_NAME MATCHES "^mpi(fort|f90)$")
+  foreach(mpi_include_dir IN LISTS MPI_Fortran_INCLUDE_DIRS)
+    if(mpi_include_dir)
+      get_filename_component(mpi_include_dir "${mpi_include_dir}" ABSOLUTE)
+      list(APPEND PETSC_MPI_INCLUDE_FLAGS "-I${mpi_include_dir}")
+    endif()
+  endforeach()
+endif()
 
+# ========================
+# PETSc (REQUIRED)
+# ========================
+add_library(HeatFlowPetsc INTERFACE)
+add_library(HeatFlow::PETSc ALIAS HeatFlowPetsc)
 
-set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${PPFLAGS}")
+find_package(PkgConfig QUIET)
 
+if(PKG_CONFIG_FOUND)
+  pkg_search_module(PETSC QUIET IMPORTED_TARGET petsc PETSc)
+endif()
 
-message(STATUS "library files: ${PREPENDED_LIB_FILES}")
-message(STATUS "source files: ${PREPENDED_SRC_FILES}")
-message(STATUS "project name: ${PROJECT_NAME}")
+if(TARGET PkgConfig::PETSC)
+  message(STATUS "Using PETSc via pkg-config")
+  target_link_libraries(HeatFlowPetsc INTERFACE PkgConfig::PETSC)
 
-# create the library
-add_library(${PROJECT_LIBRARY_NAME} STATIC ${PREPENDED_LIB_FILES})
-set_target_properties(${PROJECT_LIBRARY_NAME} PROPERTIES Fortran_MODULE_DIRECTORY ${MODULE_DIR})
-target_link_libraries(${PROJECT_LIBRARY_NAME} PUBLIC)
-target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Release>:${OPTIMFLAGS}>")
-target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Bigmem>:${MEMFLAGS}>")
-target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Parallel>:${OPTIMFLAGS}>")
-target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Parallel>:${MPFLAGS}>")
-target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Dev>:${DEVFLAGS}>")
-target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Debug>:${DEBUGFLAGS}>")
-target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Debug>:${WARNFLAGS}>")
-if(NOT CMAKE_SYSTEM_NAME STREQUAL "Darwin")
-    target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Release>:${MEMFLAGS}>")
-    target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Parallel>:${MEMFLAGS}>")
-    target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Dev>:${MEMFLAGS}>")
-    target_compile_options(${PROJECT_LIBRARY_NAME} PUBLIC "$<$<CONFIG:Debug>:${MEMFLAGS}>")
-endif()
+else()
+  if(NOT DEFINED PETSC_DIR)
+    set(PETSC_DIR $ENV{PETSC_DIR})
+  endif()
 
-add_executable(${PROJECT_NAME} ${PREPENDED_SRC_FILES})
-target_link_libraries(${PROJECT_NAME} PRIVATE ${PROJECT_LIBRARY_NAME})
-install(TARGETS ${PROJECT_NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
-set_target_properties(${PROJECT_NAME} PROPERTIES Fortran_MODULE_DIRECTORY ${MODULE_DIR})
-
-# set compile options based on different build configurations
-target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Release>:${OPTIMFLAGS}>")
-target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Bigmem>:${MEMFLAGS}>")
-target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Parallel>:${OPTIMFLAGS}>")
-target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Parallel>:${MPFLAGS}>")
-target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Dev>:${DEVFLAGS}>")
-target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Debug>:${DEBUGFLAGS}>")
-target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Debug>:${WARNFLAGS}>")
-if(NOT CMAKE_SYSTEM_NAME STREQUAL "Darwin")
-    target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Release>:${MEMFLAGS}>")
-    target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Parallel>:${MEMFLAGS}>")
-    target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Dev>:${MEMFLAGS}>")
-    target_compile_options(${PROJECT_NAME} PUBLIC "$<$<CONFIG:Debug>:${MEMFLAGS}>")
-endif()
+  if(NOT DEFINED PETSC_ARCH)
+    set(PETSC_ARCH $ENV{PETSC_ARCH})
+  endif()
 
+  if(NOT PETSC_DIR OR NOT PETSC_ARCH)
+    message(FATAL_ERROR "PETSc not found. Set PETSC_DIR and PETSC_ARCH.")
+  endif()
 
-# include the build test directory
-if(BUILD_TESTS)
-  add_subdirectory(test)
+  set(PETSC_ROOT_DIR "${PETSC_DIR}")
+  set(PETSC_INCLUDE_DIRS
+    "${PETSC_ROOT_DIR}/include"
+    "${PETSC_ROOT_DIR}/${PETSC_ARCH}/include")
+
+  find_library(PETSC_LIBRARY
+    NAMES petsc
+    HINTS "${PETSC_ROOT_DIR}/${PETSC_ARCH}/lib"
+    REQUIRED)
+
+  target_include_directories(HeatFlowPetsc INTERFACE ${PETSC_INCLUDE_DIRS})
+  target_link_libraries(HeatFlowPetsc INTERFACE ${PETSC_LIBRARY})
 endif()
 
-if ( ( CMAKE_Fortran_COMPILER MATCHES ".*gfortran.*" OR CMAKE_Fortran_COMPILER MATCHES ".*gcc.*" ) AND 
-     ( CMAKE_BUILD_TYPE MATCHES "Coverage" ) )
+set(PETSC_SOLVER_INCLUDE_FLAGS)
 
-    append_coverage_compiler_flags()
+foreach(petsc_include_dir IN LISTS PETSC_INCLUDE_DIRS)
+  if(petsc_include_dir)
+    get_filename_component(petsc_include_dir "${petsc_include_dir}" ABSOLUTE)
+    list(APPEND PETSC_SOLVER_INCLUDE_FLAGS "-I${petsc_include_dir}")
+  endif()
+endforeach()
 
+# ========================
+# PETSc solver (FIXED STALL ROOT CAUSE)
+# ========================
+add_custom_command(
+  OUTPUT ${PETSC_SOLVER_OBJECT}
+  BYPRODUCTS ${MODULE_DIR}/petsc_solver.mod
+  COMMAND ${CMAKE_COMMAND} -E make_directory ${MODULE_DIR}
+  COMMAND ${CMAKE_COMMAND} -E env
+          FC=${CMAKE_Fortran_COMPILER}
+          ${CMAKE_Fortran_COMPILER}
+          ${PPFLAGS}
+          ${PETSC_MPI_INCLUDE_FLAGS}
+          ${PETSC_SOLVER_INCLUDE_FLAGS}
+          -J${MODULE_DIR}
+          -I${MODULE_DIR}
+          -c ${PETSC_SOLVER_FILE}
+          -o ${PETSC_SOLVER_OBJECT}
+  DEPENDS ${PETSC_SOLVER_FILE}
+  VERBATIM)
+
+add_custom_target(petsc_solver_object DEPENDS ${PETSC_SOLVER_OBJECT})
+set_source_files_properties(${PETSC_SOLVER_OBJECT} PROPERTIES GENERATED TRUE EXTERNAL_OBJECT TRUE)
+
+# ========================
+# Main library
+# ========================
+add_library(${PROJECT_LIBRARY_NAME} STATIC)
+
+target_sources(${PROJECT_LIBRARY_NAME}
+  PRIVATE ${PREPENDED_LIB_FILES} ${PETSC_SOLVER_OBJECT})
+
+add_dependencies(${PROJECT_LIBRARY_NAME} petsc_solver_object)
+
+set_target_properties(${PROJECT_LIBRARY_NAME}
+  PROPERTIES Fortran_MODULE_DIRECTORY ${MODULE_DIR})
+
+target_include_directories(${PROJECT_LIBRARY_NAME}
+  PUBLIC ${MODULE_DIR})
+
+target_compile_options(${PROJECT_LIBRARY_NAME}
+  PRIVATE ${PPFLAGS})
+
+target_link_libraries(${PROJECT_LIBRARY_NAME}
+  PUBLIC HeatFlow::MPI HeatFlow::PETSc)
+
+# ========================
+# Executable
+# ========================
+add_executable(${PROJECT_NAME})
+
+target_sources(${PROJECT_NAME}
+  PRIVATE ${PREPENDED_SRC_FILES})
+
+set_target_properties(${PROJECT_NAME}
+  PROPERTIES Fortran_MODULE_DIRECTORY ${MODULE_DIR})
+
+target_link_libraries(${PROJECT_NAME}
+  PRIVATE ${PROJECT_LIBRARY_NAME})
+
+# ========================
+# Coverage
+# ========================
+if (CMAKE_BUILD_TYPE MATCHES "Coverage")
+  list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake")
+  include(CodeCoverage)
+  append_coverage_compiler_flags()
 endif()
 
+# ========================
+# Testing
+# ========================
+enable_testing()
+
 if(BUILD_TESTS)
-    # Create a directory for test data
-    add_custom_command(
-        OUTPUT ${CMAKE_BUILD_DIR}/test/test_dir_created  # Temporary marker file to indicate the directory was created
-        COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_BUILD_DIR}/test/test
-        COMMENT "Creating test directory"
-    )
-
-    # Copy test data files
-    add_custom_command(
-        OUTPUT ${CMAKE_BUILD_DIR}/test/test_data_copied  # Temporary marker file to indicate the data was copied
-        COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_LIST_DIR}/test/data ${CMAKE_BUILD_DIR}/test/test/data
-        DEPENDS ${CMAKE_CURRENT_LIST_DIR}/test/data ${CMAKE_BUILD_DIR}/test/test_dir_created  # Ensure the source data exists and the directory was created
-        COMMENT "Copying test data files"
-    )
-
-    # Create a custom target to ensure that the copy commands are executed
-    add_custom_target(copy_test_data ALL
-        DEPENDS ${CMAKE_BUILD_DIR}/test/test_data_copied  # Ensure it depends on the data copying command
-    )
-
-    # Add dependency to your main project target
-    add_dependencies(${PROJECT_NAME} copy_test_data)
-
-    # Ensure that ctest waits for copy_test_data
-    add_custom_command(TARGET copy_test_data POST_BUILD
-        COMMAND ${CMAKE_COMMAND} -E echo "Test data files copied successfully"
-    )
-
-    # Get all registered tests
-    get_property(all_tests GLOBAL PROPERTY TESTS)
-
-    # Loop through all tests and add the dependency
-    foreach(test_name IN LISTS all_tests)
-        add_dependencies(${test_name} copy_test_data)
-    endforeach()
-
-endif()
\ No newline at end of file
+  add_subdirectory(test)
+endif()
+
+# ========================
+# Install
+# ========================
+install(TARGETS ${PROJECT_NAME}
+  DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
\ No newline at end of file
diff --git a/README.md b/README.md
index 1dcc697..729cebc 100644
--- a/README.md
+++ b/README.md
@@ -7,89 +7,208 @@
 
 by Harry Mclean, Francis Huw Davies, Ned Thaddeus Taylor, and Steven Paul Hepplestone
 
-HeatFlow is Fortran-based a software package for modelling dynamical heat transport in systems using finite difference methods.
-The software is primarily designed to utilise the Cattaneo method.
-However, the Fourier method can be used instead.
+HeatFlow is a Fortran-based software package for modelling dynamical heat transport in systems using finite difference methods.
+The software is primarily designed to utilise the Cattaneo method, although the Fourier method can also be used.
 
-**IMPORTANT NOTICE: Repository Migration to GitHub**
-
-Dear users and contributors,
-
-This repository has been to be migrated from the University of Exeter GitLab to GitHub to facilitate community interaction and support.
-The latest version, updates, and collaboration now take place on this GitHub repository.
+---
 
-**GitLab Repository (Archived):** https://git.exeter.ac.uk/hepplestone/heatflow-mk2
+## IMPORTANT NOTICE: Repository Migration to GitHub
 
-## Why the Migration?
+This repository has been migrated from the University of Exeter GitLab to GitHub to facilitate community interaction and support.
+The latest version, updates, and collaboration now take place here.
 
-It was decided that this project should be migrated to allow for better community support (i.e. allowing community users to raise issues).
-All information has been ported over where possible.
-Releases prior to `HeatFlow_CattaneoPaper` have had their history modified to remove history of files over 50MB in size.
+**GitLab Repository (Archived):**
+[https://git.exeter.ac.uk/hepplestone/heatflow-mk2](https://git.exeter.ac.uk/hepplestone/heatflow-mk2)
 
-## How to Contribute on GitHub?
+### Why the Migration?
 
-Thank you for your understanding and continued support!
+The move enables better community support, including issue tracking and collaboration.
+All information has been ported where possible.
+Releases prior to `HeatFlow_CattaneoPaper` have had their history modified to remove files larger than 50 MB.
 
 ---
 
-
 ## Requirements
 
-- Fortran compiler supporting Fortran 2003 standard or later
-- fpm or CMake
+* Fortran compiler supporting Fortran 2003 or later
+* fpm or CMake
+* PETSc
+* MPI
+
+Tested with:
 
-The software bas been developed and tested using the following Fortran compilers:
-- gfortran -- gcc 13.2.0
-- gfortran -- gcc 14.1.0
+* gfortran (GCC 13.2.0)
+* gfortran (GCC 14.1.0)
+
+---
 
 ## Installation
 
-To install HeatFlow, the source must be obtained from the git repository. Use the following commands to get started:
+First obtain the source:
+
+```bash
+git clone https://github.com/ExeQuantCode/HeatFlow.git
+cd HeatFlow
 ```
- git clone https://github.com/ExeQuantCode/HeatFlow.git
- cd HeatFlow
+
+---
+
+## Building with fpm
+
+PETSc is configured at build time (not hard-coded), making the build portable across systems.
+
+### 1. Set PETSc environment variables
+
+```bash
+export PETSC_DIR=/path/to/petsc
+export PETSC_ARCH=arch-your-build
 ```
 
-### fpm
+These are standard variables provided by PETSc.
 
-To install using fpm, run the following command in the repository root directory:
+---
+
+### 2. Build
 
+```bash
+fpm build \
+  --flag "-I${PETSC_DIR}/include -I${PETSC_DIR}/${PETSC_ARCH}/include" \
+  --link-flag "-L${PETSC_DIR}/${PETSC_ARCH}/lib -Wl,-rpath,${PETSC_DIR}/${PETSC_ARCH}/lib"
 ```
-fpm build --profile=release
+
+---
+
+### 3. Run
+
+```bash
+fpm run --profile petsc -- [ALL PROGRAM OPTIONS]
 ```
 
-To execute the code, use
+---
+
+### Alternative: pkg-config
+
+If PETSc provides pkg-config:
 
+```bash
+fpm build \
+  --flag "$(pkg-config --cflags petsc)" \
+  --link-flag "$(pkg-config --libs petsc)"
 ```
-fpm run HeatFlow --profile release -- [ALL PROGRAM OPTIONS]
+
+---
+
+### Notes (fpm)
+
+* No PETSc paths are stored in `fpm.toml`
+* All system-specific configuration is provided at build time
+* Works well on clusters, macOS, and CI environments
+
+---
+
+## Building with CMake
+
+CMake provides automatic detection of MPI and PETSc.
+
+### 1. Configure PETSc
+
+If PETSc is not installed system-wide:
+
+```bash
+export PETSC_DIR=/path/to/petsc
+export PETSC_ARCH=arch-your-build
 ```
 
-### cmake
+---
 
-For cmake installation, start within the repository root directory, run the following commands:
+### 2. Configure and build
 
+```bash
+cmake -B build -S . -DCMAKE_BUILD_TYPE=Release
+cmake --build build
 ```
-mkdir build
-cd build
-cmake [-DCMAKE_BUILD_TYPE=Release] ..
-make install
+
+Or explicitly:
+
+```bash
+cmake -B build -S . \
+  -DCMAKE_BUILD_TYPE=Release \
+  -DPETSC_DIR=$PETSC_DIR \
+  -DPETSC_ARCH=$PETSC_ARCH
 ```
 
-This will build and install the executable in the following directory:
+---
+
+### 3. Install
+
+```bash
+cmake --install build
 ```
+
+This installs the executable to:
+
+```bash
 ${HOME}/.local/HeatFlow/bin/HeatFlow
 ```
 
-This executable can now be called to run the HeatFlow software package and simulate heat transport.
-If the `${HOME}/.local/HeatFlow/bin` is added to your `PATH` environment variable, then the program can be called as a terminal command.
-This can be done with the following command (works on a per-terminal basis, if you want to update it for all, include this in your source shell file):
+Add to your `PATH` if desired:
 
-```
+```bash
 export PATH="${PATH}:${HOME}/.local/HeatFlow/bin"
 ```
 
-To execute the program, use the following command:
+---
 
-```
+### 4. Run
+
+```bash
 HeatFlow
 ```
+
+Optional directory flags:
+
+```bash
+HeatFlow --input-directory /path/to/run/inputs --output-directory /path/to/run/outputs
+```
+
+Compatibility alias for the legacy run-layout:
+
+```bash
+HeatFlow --directory /path/to/run
+```
+
+The legacy alias maps to `/path/to/run/inputs`, `/path/to/run/outputs`, and `/path/to/run/restart`.
+
+---
+
+### Using pkg-config (recommended)
+
+If PETSc supports pkg-config, it will be detected automatically:
+
+```bash
+cmake -B build -S .
+cmake --build build
+```
+
+---
+
+### Notes (CMake)
+
+* PETSc detection order:
+
+  1. pkg-config
+  2. `PETSC_DIR` / `PETSC_ARCH`
+  3. Local `../petsc` fallback
+
+* MPI is detected automatically
+
+* Configuration fails with a clear error if PETSc is not found
+
+---
+
+## Summary
+
+* **fpm**: lightweight and flexible, ideal for development
+* **CMake**: robust and portable, suited for deployment
+
+In both cases, PETSc is configured at build time rather than hard-coded, ensuring portability across systems.
diff --git a/build_fpm_with_petsc.sh b/build_fpm_with_petsc.sh
new file mode 100755
index 0000000..d201feb
--- /dev/null
+++ b/build_fpm_with_petsc.sh
@@ -0,0 +1,8 @@
+#!/usr/bin/env bash
+
+: "${PETSC_DIR:?Need to set PETSC_DIR}"
+: "${PETSC_ARCH:?Need to set PETSC_ARCH}"
+
+fpm build \
+  --flag "-I${PETSC_DIR}/include -I${PETSC_DIR}/${PETSC_ARCH}/include" \
+  --link-flag "-L${PETSC_DIR}/${PETSC_ARCH}/lib -Wl,-rpath,${PETSC_DIR}/${PETSC_ARCH}/lib"
\ No newline at end of file
diff --git a/fpm.toml b/fpm.toml
index fcc5f3f..0b50c81 100644
--- a/fpm.toml
+++ b/fpm.toml
@@ -15,8 +15,6 @@ 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"]
diff --git a/src/heatflow.f90 b/src/heatflow.f90
index cc16531..a0c306b 100644
--- a/src/heatflow.f90
+++ b/src/heatflow.f90
@@ -22,9 +22,10 @@ program HEATFLOW_V0_3
   
   use constants, only: real12, int12
   use constructions, only: heatblock
-  use output, only: data_write, final_print
-  use inputs, only: read_all_files, iverb, ntime, LPercentage
-  use inputs, only: IVERB
+   use output, only: data_write, final_print
+   use inputs, only: read_all_files, iverb, ntime, LPercentage
+   use inputs, only: IVERB, input_directory, output_directory, restart_directory
+   use inputs, only: set_io_directories, join_path
   use evolution, only: simulate
   use setup, only: set_global_variables
   use INITIAL, only: initial_evolve
@@ -35,7 +36,7 @@ program HEATFLOW_V0_3
    integer(int12) :: itime
 
    !-------------------------------------------------------------!
-   ! Initialize PETSc FIRST (before any other operations)       !
+   ! Initialize PETSc FIRST (before any other operations)        !
    !-------------------------------------------------------------!
    CALL petsc_init()
    !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
@@ -47,19 +48,26 @@ program HEATFLOW_V0_3
    !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
 
    ! give feedback to user that code has begun
-   if (petsc_is_root()) write(*,*) 'Setup initialising'
-   
+   if(petsc_is_root()) write(*,*) 'Setup initialising'
+
+   !-------------------------------------------------------------!
+   ! handle command line arguments on master/root node only      !
+   !-------------------------------------------------------------!
+   call handle_command_line_arguments()
+
    !-------------------------------------------------------------!
    ! Read parameters from input file and set global variables ...!
    ! ... and arrays                                              !
    !-------------------------------------------------------------!
-   CALL read_all_files()                                         
-   
-   CALL cpu_time(cpustart2)                                      
-   CALL set_global_variables() 
-   CALL cpu_time(cpuend)
-   if (IVERB.ge.1) write(*,'(A,F12.6)') &
-   ' time to complete set_global_variables=', cpuend-cpustart2   
+   if(petsc_is_root())then
+      CALL read_all_files()                                         
+      
+      CALL cpu_time(cpustart2)                                      
+      CALL set_global_variables() 
+      CALL cpu_time(cpuend)
+      if (IVERB.ge.1) write(*,'(A,F12.6)') &
+      ' time to complete set_global_variables=', cpuend-cpustart2   
+   end if
 
  
    !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
@@ -111,6 +119,101 @@ program HEATFLOW_V0_3
    ! give feedback to user that code has ended
    if (petsc_is_root()) write(*,*) 'all done'
 
+contains
+
+   !-------------------------------------------------------------!
+   ! Handle command line arguments for directory specification   !
+   !-------------------------------------------------------------!
+   subroutine handle_command_line_arguments()
+     implicit none
+     integer :: nargs, i
+     character(len=1024) :: arg, directory, cli_input_directory, cli_output_directory
+     logical :: dir_exists
+     character(len=*), parameter :: directory_flag = '--directory'
+     character(len=*), parameter :: directory_prefix = '--directory='
+     character(len=*), parameter :: input_directory_flag = '--input-directory'
+     character(len=*), parameter :: input_directory_prefix = '--input-directory='
+     character(len=*), parameter :: output_directory_flag = '--output-directory'
+     character(len=*), parameter :: output_directory_prefix = '--output-directory='
+      
+     nargs = command_argument_count()
+     directory = ''
+     cli_input_directory = ''
+     cli_output_directory = ''
+      
+     if (petsc_is_root()) write(*,*) 'Number of command line arguments: ', nargs
+     i = 1
+     do while (i .le. nargs)
+        call get_command_argument(i, arg)
+        if (petsc_is_root()) write(*,*) 'Received command line argument: ', trim(arg)
+        if (trim(arg) .eq. directory_flag) then
+           call require_argument_value(i, nargs, directory_flag, directory)
+           i = i + 2
+           cycle
+        else if (index(trim(arg), directory_prefix) .eq. 1) then
+           directory = trim(arg(len(directory_prefix) + 1:))
+        else if (trim(arg) .eq. input_directory_flag) then
+           call require_argument_value(i, nargs, input_directory_flag, cli_input_directory)
+           i = i + 2
+           cycle
+        else if (index(trim(arg), input_directory_prefix) .eq. 1) then
+           cli_input_directory = trim(arg(len(input_directory_prefix) + 1:))
+        else if (trim(arg) .eq. output_directory_flag) then
+           call require_argument_value(i, nargs, output_directory_flag, cli_output_directory)
+           i = i + 2
+           cycle
+        else if (index(trim(arg), output_directory_prefix) .eq. 1) then
+           cli_output_directory = trim(arg(len(output_directory_prefix) + 1:))
+        end if
+        i = i + 1
+     end do
+
+     if(len_trim(directory) .gt. 0) then
+        inquire(file=trim(directory)//'/.' , exist=dir_exists)
+        if(.not. dir_exists) then
+           if (petsc_is_root()) write(*,*) 'Error: Directory does not exist: ', trim(directory)
+           call exit(1)
+        end if
+        call set_io_directories(input_dir=join_path(directory, 'inputs'), &
+             output_dir=join_path(directory, 'outputs'), restart_dir=join_path(directory, 'restart'))
+     end if
+
+     if(len_trim(cli_input_directory) .gt. 0) then
+        inquire(file=trim(cli_input_directory)//'/.' , exist=dir_exists)
+        if(.not. dir_exists) then
+           if (petsc_is_root()) write(*,*) 'Error: Input directory does not exist: ', trim(cli_input_directory)
+           call exit(1)
+        end if
+        call set_io_directories(input_dir=cli_input_directory)
+     end if
+
+     if(len_trim(cli_output_directory) .gt. 0) then
+        call set_io_directories(output_dir=cli_output_directory)
+     end if
+
+     if (petsc_is_root()) then
+        write(*,*) 'Input directory: ', trim(input_directory)
+        write(*,*) 'Output directory: ', trim(output_directory)
+        write(*,*) 'Restart directory: ', trim(restart_directory)
+     end if
+
+   end subroutine handle_command_line_arguments
+   !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
+
+     subroutine require_argument_value(index, count, flag_name, value)
+       implicit none
+       integer, intent(in) :: index, count
+       character(len=*), intent(in) :: flag_name
+       character(len=*), intent(out) :: value
+
+       if (index .eq. count) then
+          if (petsc_is_root()) write(*,*) 'Error: Missing value for ', trim(flag_name)
+          call exit(1)
+       end if
+
+       call get_command_argument(index + 1, value)
+     end subroutine require_argument_value
+
 end program HEATFLOW_V0_3
 
 
diff --git a/src/heatflow/mod_init_evolve.f90 b/src/heatflow/mod_init_evolve.f90
index b23144f..ff038d5 100644
--- a/src/heatflow/mod_init_evolve.f90
+++ b/src/heatflow/mod_init_evolve.f90
@@ -21,6 +21,7 @@
 module initial
    use constants, only: real12, int12, fields
    use inputs, only: NA, nx, ny, nz, InputTempDis, FullRestart, T_System
+   use inputs, only: restart_directory, join_path
    use globe_data, only: Temp_p,Temp_pp
    implicit none
 
@@ -38,10 +39,10 @@ subroutine initial_evolve()
       !------------------------------------------------------------------------------
 
       if (FullRestart) then
-         CALL read_temp_file('./restart/TempDis.dat',Temp_p)
-         CALL read_temp_file('./restart/TempDisTPD.dat',Temp_pp)
+         CALL read_temp_file(join_path(restart_directory, 'TempDis.dat'),Temp_p)
+         CALL read_temp_file(join_path(restart_directory, 'TempDisTPD.dat'),Temp_pp)
       else if (InputTempDis) then
-         CALL read_temp_file('./restart/TempDis.dat',Temp_p)
+         CALL read_temp_file(join_path(restart_directory, 'TempDis.dat'),Temp_p)
          Temp_pp = Temp_p
       else
          Temp_p = T_System
diff --git a/src/heatflow/mod_inputs.f90 b/src/heatflow/mod_inputs.f90
index 7f1fd9b..edc5d7b 100644
--- a/src/heatflow/mod_inputs.f90
+++ b/src/heatflow/mod_inputs.f90
@@ -97,6 +97,9 @@ module inputs
 
   ! Name of simiulation run
   character(1024) :: RunName
+  character(1024) :: input_directory = './inputs'
+  character(1024) :: output_directory = './outputs'
+  character(1024) :: restart_directory = './restart'
   character(12)::Periodic
   ! Essentially it is the system that is being simulated
   type(heatblock), dimension(:,:,:), allocatable :: grid 
@@ -111,13 +114,94 @@ module inputs
 
 
 contains
+!!!#################################################################################################
+!!! configure the directories used for inputs, outputs, and restart files
+!!!#################################################################################################
+  subroutine set_io_directories(input_dir, output_dir, restart_dir)
+   implicit none
+   character(len=*), intent(in), optional :: input_dir, output_dir, restart_dir
+
+   if (present(input_dir)) then
+     if (len_trim(input_dir) .gt. 0) then
+       input_directory = normalize_directory(input_dir)
+       call derive_restart_directory_from_input()
+     end if
+   end if
+
+   if (present(output_dir)) then
+     if (len_trim(output_dir) .gt. 0) output_directory = normalize_directory(output_dir)
+   end if
+
+   if (present(restart_dir)) then
+     if (len_trim(restart_dir) .gt. 0) restart_directory = normalize_directory(restart_dir)
+   end if
+  end subroutine set_io_directories
+!!!#################################################################################################
+
+!!!#################################################################################################
+!!! build a filesystem path from a directory and file name
+!!!#################################################################################################
+  function join_path(directory, name) result(path)
+   implicit none
+   character(len=*), intent(in) :: directory, name
+   character(len=1024) :: path
+
+   if (len_trim(directory) .eq. 0) then
+     path = trim(name)
+   else if (len_trim(name) .eq. 0) then
+     path = trim(directory)
+   else
+     path = trim(normalize_directory(directory)) // '/' // trim(name)
+   end if
+  end function join_path
+!!!#################################################################################################
+
+!!!#################################################################################################
+!!! remove trailing slashes so path joins stay consistent
+!!!#################################################################################################
+  function normalize_directory(dirname) result(normalized)
+   implicit none
+   character(len=*), intent(in) :: dirname
+   character(len=1024) :: normalized
+   integer :: last_char
+
+   normalized = trim(adjustl(dirname))
+   last_char = len_trim(normalized)
+
+   do while (last_char .gt. 1 .and. normalized(last_char:last_char) .eq. '/')
+     last_char = last_char - 1
+   end do
+
+   normalized = normalized(:last_char)
+  end function normalize_directory
+!!!#################################################################################################
+
+!!!#################################################################################################
+!!! keep restart files beside an inputs directory when possible
+!!!#################################################################################################
+  subroutine derive_restart_directory_from_input()
+   implicit none
+   integer :: dir_length
+
+   dir_length = len_trim(input_directory)
+
+   if (trim(input_directory) .eq. 'inputs') then
+     restart_directory = 'restart'
+   else if (dir_length .gt. 7) then
+     if (input_directory(dir_length-6:dir_length) .eq. '/inputs') then
+       restart_directory = trim(input_directory(:dir_length-7)) // '/restart'
+     end if
+   end if
+  end subroutine derive_restart_directory_from_input
+!!!#################################################################################################
+
 !!!#################################################################################################
 !!! read_all_files will call the routines to read each input file
 !!!#################################################################################################
   subroutine read_all_files()
     implicit none
     integer :: unit, reason ! file unit and reason
-    character(64) :: param_infile, mat_infile, mesh_infile ! file names
+   character(1024) :: param_infile, mat_infile, mesh_infile ! file names
     logical :: file_exists ! check if file exists
 
 
@@ -125,7 +209,7 @@ subroutine read_all_files()
     ! get data from param.in
     !-----------------------------------------------
     ! name infile
-    param_infile = "./inputs/param.in" ! file name
+    param_infile = join_path(input_directory, 'param.in') ! file name
 
     ! check if file is there
     inquire(file=param_infile, exist=file_exists) ! check if file exists
@@ -150,7 +234,7 @@ subroutine read_all_files()
     !-----------------------------------------------
     ! name infile
 
-    mat_infile = "./inputs/mat.in" ! file name
+    mat_infile = join_path(input_directory, 'mat.in') ! file name
 
     ! check if file is there
     inquire(file=mat_infile, exist=file_exists)
@@ -173,7 +257,7 @@ subroutine read_all_files()
     !-----------------------------------------------
     ! name infile
 
-    mesh_infile = "./inputs/system.in" ! file name
+    mesh_infile = join_path(input_directory, 'system.in') ! file name
 
     ! check if file is there
     inquire(file=mesh_infile, exist=file_exists)
diff --git a/src/heatflow/mod_output.f90 b/src/heatflow/mod_output.f90
index c518761..c7c6de4 100644
--- a/src/heatflow/mod_output.f90
+++ b/src/heatflow/mod_output.f90
@@ -42,10 +42,11 @@
 !!! Author: Harry Mclean, Frank Davies, Steven Hepplestone
 !!!#################################################################################################
 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, CompressedOutput
-  use inputs, only: start_ix, end_ix, start_iy, end_iy, start_iz, end_iz
+   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, CompressedOutput
+   use inputs, only: start_ix, end_ix, start_iy, end_iy, start_iz, end_iz
+   use inputs, only: output_directory, join_path
   use globe_data, only: Temp_p,Temp_pp, heat, heated_volume, logname
   implicit none
   
@@ -60,8 +61,9 @@ subroutine data_write(itime)
     
     logunit = 20
     file_prefix = 'Temperture_'
-    outdir='./outputs/'
+   outdir = output_directory
     file_extension = '.out'
+   call ensure_directory(outdir)
     
     !---------------------------------------
     !  make a 3d array
@@ -84,8 +86,8 @@ subroutine data_write(itime)
           !---------------------------------------
           ! open test output files                
           !---------------------------------------
-          open(unit=33,file='./outputs/Power.txt')
-          open(unit=30, file='./outputs/Test.txt')
+          open(unit=33,file=join_path(outdir, 'Power.txt'))
+          open(unit=30, file=join_path(outdir, 'Test.txt'))
           !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        else
           !---------------------------------------
@@ -173,6 +175,17 @@ subroutine data_write(itime)
  end subroutine data_write
 !!!########################################################################
 
+!!!########################################################################
+ subroutine ensure_directory(dirname)
+    character(len=*), intent(in) :: dirname
+    logical :: exists
+
+   inquire(file=trim(dirname), exist=exists)
+   if (.not. exists) call execute_command_line('mkdir -p "' // trim(dirname) // '"')
+ end subroutine ensure_directory
+
+!!!########################################################################
+
 
  
 !!!########################################################################
@@ -195,7 +208,7 @@ subroutine PlotdeltaT(itime)
       end do
    end do
    if (itime .eq. 1) then
-      open(unit=31,file='./outputs/DTemperature.txt')
+      open(unit=31,file=join_path(output_directory, 'DTemperature.txt'))
    end if
    
 
@@ -210,17 +223,17 @@ subroutine last_log(logname,outdir)
    character(len=1024), intent(in) :: outdir
    logical :: flag
    integer(int12) :: i
+    character(len=1024) :: basename
    
    i = 0
    flag=.true.
    do while (flag)
       if (CompressedOutput) then
-         write(logname, '(A,A,I2.2,A)') trim(adjustl(outdir)) // 'output_' // &
-              trim(adjustl(RunName)),'_',  i, '.bin'
+         write(basename, '(A,I2.2,A)') 'output_' // trim(adjustl(RunName)) // '_', i, '.bin'
       else
-         write(logname, '(A,A,I2.2)') trim(adjustl(outdir)) // 'output_' // &
-              trim(adjustl(RunName)),'_',  i
+         write(basename, '(A,I2.2)') 'output_' // trim(adjustl(RunName)) // '_', i
       endif
+      logname = join_path(outdir, trim(adjustl(basename)))
       inquire(file=logname, exist=flag)
       i = i+1
    end do
@@ -234,6 +247,8 @@ subroutine final_print()
    real(real12) :: TotalPower, totaltime, vol
    integer(int12) :: unit
    character(len=64) :: form
+
+   call ensure_directory(output_directory)
    
    if (IVERB .gt. 3) then 
       TotalPower=heat
@@ -248,14 +263,14 @@ subroutine final_print()
    
    end if 
 
-    open(newunit=unit,file='./outputs/TempDis.dat')
+   open(newunit=unit,file=join_path(output_directory, 'TempDis.dat'))
     write(form,'(A,I0,A)') '(',fields,'(ES16.8,1X))'
     write(unit,form) Temp_p(:)
     write(unit,*)
     write(unit,*)
     close(unit)
     
-    open(newunit=unit,file='./outputs/TempDisTPD.dat')
+   open(newunit=unit,file=join_path(output_directory, 'TempDisTPD.dat'))
     write(form,'(A,I0,A)') '(',fields,'(ES16.8,1X))'
     write(unit,form) Temp_pp(:)
     write(unit,*)
diff --git a/src/heatflow/mod_petsc_solver.f90 b/src/petsc_solver/mod_petsc_solver.f90
similarity index 94%
rename from src/heatflow/mod_petsc_solver.f90
rename to src/petsc_solver/mod_petsc_solver.f90
index 65bbf45..310bc65 100644
--- a/src/heatflow/mod_petsc_solver.f90
+++ b/src/petsc_solver/mod_petsc_solver.f90
@@ -16,10 +16,10 @@ module petsc_solver
   ! ====================================
 
   ! 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
+  Mat, save :: A_saved
+  Vec, save :: bb_saved
+  Vec, save :: xx_saved
+  KSP, save :: ksp_saved
   logical, save :: initialized = .false.
   integer, save :: n_saved = 0
   integer, save :: comm_rank_saved = 0
@@ -38,6 +38,10 @@ module petsc_solver
   subroutine petsc_init()
     integer :: ierr
     call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
+    PetscObjectNullify(A_saved)
+    PetscObjectNullify(bb_saved)
+    PetscObjectNullify(xx_saved)
+    PetscObjectNullify(ksp_saved)
     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
@@ -59,10 +63,22 @@ 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 (.not. PetscObjectIsNull(A_saved)) then
+      call MatDestroy(A_saved, ierr)
+      PetscObjectNullify(A_saved)
+    end if
+    if (.not. PetscObjectIsNull(bb_saved)) then
+      call VecDestroy(bb_saved, ierr)
+      PetscObjectNullify(bb_saved)
+    end if
+    if (.not. PetscObjectIsNull(xx_saved)) then
+      call VecDestroy(xx_saved, ierr)
+      PetscObjectNullify(xx_saved)
+    end if
+    if (.not. PetscObjectIsNull(ksp_saved)) then
+      call KSPDestroy(ksp_saved, ierr)
+      PetscObjectNullify(ksp_saved)
+    end if
     if (allocated(ia_saved)) deallocate(ia_saved)
     if (allocated(ja_saved)) deallocate(ja_saved)
     if (allocated(local_indices_saved)) deallocate(local_indices_saved)
@@ -74,10 +90,6 @@ subroutine petsc_cleanup()
     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
diff --git a/test/test_mod_material.f90 b/test/test_mod_material.f90
index bc4820e..2ebc997 100644
--- a/test/test_mod_material.f90
+++ b/test/test_mod_material.f90
@@ -6,6 +6,7 @@ program test_mod_material
 
     integer(int12) :: imaterial_type
     real(real12) :: kappa, kappa3D, h_conv, heat_capacity, rho, sound_speed, tau, em
+    real(real12), dimension(3) :: vel
     logical :: test_passed
     integer :: unit, reason, i
     logical :: file_exists
@@ -16,7 +17,8 @@ program test_mod_material
     call getcwd(cwdstring)
     ! Test case 1: Silicon (140)
     imaterial_type = 140
-    call material(imaterial_type, kappa, kappa3D, h_conv, heat_capacity, rho, sound_speed, tau, em)
+    vel = 0.0_real12
+    call material(imaterial_type, kappa, kappa3D, h_conv, heat_capacity, rho, sound_speed, tau, em, vel)
     
     test_passed = abs(kappa - 130.0_real12) < 1e-6 .and. &
                                 abs(kappa3D - 130.0_real12) < 1e-6 .and. &
@@ -37,7 +39,8 @@ program test_mod_material
     ! Test case 2: All materials
     do i = 1, size(test_mats)
         imaterial_type = test_mats(i)
-        call material(imaterial_type, kappa, kappa3D, h_conv, heat_capacity, rho, sound_speed, tau, em)
+        vel = 0.0_real12
+        call material(imaterial_type, kappa, kappa3D, h_conv, heat_capacity, rho, sound_speed, tau, em, vel)
         ! Program should exit with error message
     end do
 
@@ -70,5 +73,6 @@ program test_mod_material
         ! Test case 4: Invalid material type
     imaterial_type = -1
     write(*,*) "Testing invalid material type (-1)..."
-    call material(imaterial_type, kappa, kappa3D, h_conv, heat_capacity, rho, sound_speed, tau, em)
+    vel = 0.0_real12
+    call material(imaterial_type, kappa, kappa3D, h_conv, heat_capacity, rho, sound_speed, tau, em, vel)
 end program test_mod_material
\ No newline at end of file