diff --git a/.gitignore b/.gitignore index e6500e5a1..24f757c5c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,36 +1,37 @@ -## C++ ######################################################################### - -# Prerequisites -*.d - -# Compiled Object files -*.slo -*.lo -*.o -*.obj - -# Precompiled Headers -*.gch -*.pch - -# Compiled Dynamic libraries -*.so -*.dylib -*.dll - -# Fortran module files -*.mod -*.smod - -# Compiled Static libraries -*.lai -*.la -*.a -*.lib - -# Executables -*.exe -*.out -*.app - -*/build/* \ No newline at end of file +## C++ ######################################################################### + +# Prerequisites +*.d + +# Compiled Object files +*.slo +*.lo +*.o +*.obj + +# Precompiled Headers +*.gch +*.pch + +# Compiled Dynamic libraries +*.so +*.dylib +*.dll + +# Fortran module files +*.mod +*.smod + +# Compiled Static libraries +*.lai +*.la +*.a +*.lib + +# Executables +*.exe +*.out +*.app + +build/ +.vscode \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 6d3567742..ac2447187 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,10 @@ -cmake_minimum_required(VERSION 3.12.0) -project(project_title LANGUAGES CXX C) - -include(../common/cmake-common.cmake) - -add_executable(lab-01 src/lab-01.cpp src/Poisson1D.cpp) -deal_ii_setup_target(lab-01) +cmake_minimum_required(VERSION 3.12.0) +project(project_title LANGUAGES CXX C) + +include(common/cmake-common.cmake) + +add_executable(H_Heat src/homogeneous/main.cpp src/homogeneous/H_Heat.cpp src/timing/Timing.cpp) +deal_ii_setup_target(H_Heat) + +add_executable(STA_Heat src/adaptive/main.cpp src/adaptive/STA_Heat.cpp src/homogeneous/H_Heat.cpp src/timing/Timing.cpp) +deal_ii_setup_target(STA_Heat) diff --git a/README.md b/README.md index b2f00611b..2f06cabcb 100644 --- a/README.md +++ b/README.md @@ -1,23 +1,78 @@ -### Organizing the source code -Please place all your sources into the `src` folder. - -Binary files must not be uploaded to the repository (including executables). - -Mesh files should not be uploaded to the repository. If applicable, upload `gmsh` scripts with suitable instructions to generate the meshes (and ideally a Makefile that runs those instructions). If not applicable, consider uploading the meshes to a different file sharing service, and providing a download link as part of the building and running instructions. - -### Compiling -To build the executable, make sure you have loaded the needed modules with -```bash -$ module load gcc-glibc dealii -``` -Then run the following commands: -```bash -$ mkdir build -$ cd build -$ cmake .. -$ make -``` -The executable will be created into `build`, and can be executed through -```bash -$ ./executable-name -``` \ No newline at end of file +# Space-Time Adaptivity + +## Project Information + +- **Course:** Numerical Methods for Partial Differential Equations +- **Authors:** Giacomo Maglio, Francesca Morciano, Giuseppe Hares +- **Advisor:** Alfio Maria Quarteroni +- **Co-Advisor:** Michele Bucelli +- **Academic year:** 2025/2026 +- **Institution:** Politecnico di Milano + +This project contains two deal.II-based solvers for the heat equation: + +- `H_Heat`: homogeneous reference solver +- `STA_Heat`: adaptive space-time solver + +The code is organized under `src/` and builds with CMake. + +## Repository layout + +- `src/homogeneous/`: reference solver implementation +- `src/adaptive/`: adaptive solver implementation +- `src/timing/`: profiling utilities + +## Requirements + +You need a working C++ toolchain and deal.II with MPI support. + +Environment setup: + +```bash +$ apptainer shell *.sif +$ source /u/sw/etc/bash.bashrc +$ module load gcc-glibc dealii +``` + +## Build + +From the repository root: + +```bash +$ mkdir build +$ cd build +$ cmake .. +$ make +``` + +This creates the executables inside `build/`. + +## Run + +Available executables: + +- `build/H_Heat` +- `build/STA_Heat` + +Run them from the build directory, for example: + +```bash +$ cd build +$ mpirun -n ./H_Heat # to run the homogenous solver +$ mpirun -n ./STA_Heat # to run the adaptive solver +``` + +Both executables currently use hardcoded parameters in their respective `main.cpp` files. +To change time step, final time, or forcing term, edit: + +- `src/homogeneous/main.cpp` +- `src/adaptive/main.cpp` + +## Output + +The solvers write output files in the build directory. + +## Notes + +- The adaptive executable contains an optional baseline-comparison path behind the `COMPARE_WITH_BASE` macro. +- If you want to compare adaptive and homogeneous results, enable that macro in `src/adaptive/main.cpp` and ensure the baseline mesh is available. diff --git a/common/cmake-common.cmake b/common/cmake-common.cmake new file mode 100644 index 000000000..8e213c37a --- /dev/null +++ b/common/cmake-common.cmake @@ -0,0 +1,34 @@ +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED "ON") +set(CMAKE_CXX_FLAGS "-O3 -march=native -ffast-math") + +# Set default build type to Release. +if(NOT CMAKE_BUILD_TYPE OR "${CMAKE_BUILD_TYPE}" STREQUAL "") + set(CMAKE_BUILD_TYPE "Release" CACHE STRING "" FORCE) +endif() +message(STATUS) +message(STATUS "Build type: ${CMAKE_BUILD_TYPE}") +message(STATUS) +if("${CMAKE_BUILD_TYPE}" STREQUAL "Debug") + add_definitions(-DBUILD_TYPE_DEBUG) +endif() + +# Locate MPI compiler. +find_package(MPI REQUIRED) +set(CMAKE_CXX_COMPILER "${MPI_CXX_COMPILER}") + +# Locate Boost. +find_package(Boost 1.72.0 REQUIRED + COMPONENTS filesystem iostreams serialization + HINTS ${BOOST_DIR} $ENV{BOOST_DIR} $ENV{mkBoostPrefix}) +message(STATUS "Using the Boost-${Boost_VERSION} configuration found at ${Boost_DIR}") +message(STATUS) +include_directories(${Boost_INCLUDE_DIRS}) + +# Locate deal.II and initialize its variables. +find_package(deal.II 9.3.1 REQUIRED + HINTS ${DEAL_II_DIR} $ENV{DEAL_II_DIR} $ENV{mkDealiiPrefix}) +deal_ii_initialize_cached_variables() + +# Add useful compiler flags. +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wfloat-conversion -Wmissing-braces -Wnon-virtual-dtor") \ No newline at end of file diff --git a/experiments_results/baseline_comparison.txt b/experiments_results/baseline_comparison.txt new file mode 100644 index 000000000..6c1812313 --- /dev/null +++ b/experiments_results/baseline_comparison.txt @@ -0,0 +1,110 @@ +======================================================================= + Space and Time Adaptivity +======================================================================= + +Baseline Results +Number of steps: 40 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +update-ghost 40 1.17e-04 1.92e-04 4.22e-04 7.69e-03 +solve 40 6.74e-02 8.87e-02 1.37e-01 3.55e+00 +output 40 7.12e-02 8.70e-02 1.39e-01 3.48e+00 +assemble 40 1.82e-01 2.13e-01 2.61e-01 8.50e+00 +setup 1 3.13e-01 3.13e-01 3.13e-01 3.13e-01 +--------------------------------------------------------------- +Number of DoFs: 35937 + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 12 1.96e-01 3.56e-01 8.15e-01 4.28e+00 +time-adap 28 1.06e-01 1.41e-01 2.12e-01 3.95e+00 +solve 42 2.27e-01 3.12e-01 4.82e-01 1.31e+01 +output 28 1.06e-01 1.40e-01 2.11e-01 3.93e+00 +assemble 42 1.87e-01 3.26e-01 5.21e-01 1.37e+01 +update-ghost 73 1.10e-05 3.20e-04 1.20e-03 2.34e-02 +setup 1 1.79e-01 1.79e-01 1.79e-01 1.79e-01 +--------------------------------------------------------------- +Max number of DoFs: 82891 + +Computing L2 against baseline +L2_error = 1.33e-03 + +======================================================================= + Only Space Adaptivity +======================================================================= + +Baseline Results +Number of steps: 40 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +update-ghost 40 1.19e-04 1.92e-04 5.98e-04 7.67e-03 +solve 40 7.31e-02 8.91e-02 1.33e-01 3.56e+00 +output 40 7.35e-02 8.89e-02 1.15e-01 3.56e+00 +assemble 40 1.84e-01 2.14e-01 2.98e-01 8.55e+00 +setup 1 3.38e-01 3.38e-01 3.38e-01 3.38e-01 +--------------------------------------------------------------- +Number of DoFs: 35937 + +Adaptivity Results +Number of steps: 20 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 5 2.21e-01 5.87e-01 1.04e+00 2.93e+00 +time-adap 20 8.22e-02 1.60e-01 2.64e-01 3.20e+00 +solve 20 1.39e-01 3.44e-01 5.44e-01 6.89e+00 +output 20 8.19e-02 1.59e-01 2.62e-01 3.17e+00 +assemble 20 2.09e-01 3.70e-01 5.98e-01 7.40e+00 +update-ghost 43 1.39e-05 3.86e-04 1.42e-03 1.66e-02 +setup 1 2.55e-01 2.55e-01 2.55e-01 2.55e-01 +--------------------------------------------------------------- +Max number of DoFs: 99978 + +Computing L2 against baseline +L2_error = 6.56e-05 + + +======================================================================= + Only Time Adaptivity +======================================================================= + +Baseline Results +Number of steps: 40 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +update-ghost 40 1.09e-04 2.17e-04 5.24e-04 8.66e-03 +solve 40 6.77e-02 8.45e-02 1.17e-01 3.38e+00 +output 40 6.84e-02 8.88e-02 1.41e-01 3.55e+00 +assemble 40 1.83e-01 2.16e-01 2.72e-01 8.63e+00 +setup 1 1.80e-01 1.80e-01 1.80e-01 1.80e-01 +--------------------------------------------------------------- +Number of DoFs: 35937 + +Adaptivity Results +Number of steps: 40 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +time-adap 27 7.31e-02 8.56e-02 1.12e-01 2.31e+00 +solve 40 1.28e-01 1.92e-01 2.57e-01 7.69e+00 +output 27 7.27e-02 8.51e-02 1.12e-01 2.30e+00 +assemble 40 1.72e-01 2.07e-01 2.49e-01 8.27e+00 +update-ghost 67 5.00e-05 2.39e-04 6.81e-04 1.60e-02 +setup 1 2.26e-01 2.26e-01 2.26e-01 2.26e-01 +--------------------------------------------------------------- +Max number of DoFs: 35937 + +Computing L2 against baseline +L2_error = 1.35e-03 \ No newline at end of file diff --git a/experiments_results/scalability_comparison.txt b/experiments_results/scalability_comparison.txt new file mode 100644 index 000000000..35a63b166 --- /dev/null +++ b/experiments_results/scalability_comparison.txt @@ -0,0 +1,376 @@ +======================================================================= + Number of processes = 1 +======================================================================= + +~ RUN 1 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 12 2.32e-01 4.82e-01 9.60e-01 5.78e+00 +time-adap 28 8.43e-02 1.47e-01 2.38e-01 4.12e+00 +solve 42 1.74e-01 2.79e-01 4.09e-01 1.17e+01 +output 28 8.39e-02 1.46e-01 2.36e-01 4.10e+00 +assemble 42 1.83e-01 3.15e-01 5.50e-01 1.32e+01 +update-ghost 73 1.23e-05 3.10e-04 1.19e-03 2.27e-02 +setup 1 2.05e-01 2.05e-01 2.05e-01 2.05e-01 +--------------------------------------------------------------- +Max number of DoFs: 82891 + +~ RUN 2 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 12 2.85e-01 5.05e-01 1.02e+00 6.06e+00 +time-adap 28 8.92e-02 1.65e-01 2.95e-01 4.62e+00 +solve 42 2.01e-01 3.23e-01 4.98e-01 1.36e+01 +output 28 8.90e-02 1.64e-01 2.93e-01 4.60e+00 +assemble 42 2.00e-01 3.49e-01 5.03e-01 1.47e+01 +update-ghost 73 1.07e-05 3.17e-04 1.39e-03 2.32e-02 +setup 1 2.40e-01 2.40e-01 2.40e-01 2.40e-01 +--------------------------------------------------------------- +Max number of DoFs: 82891 + +~ RUN 3 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 12 2.62e-01 5.86e-01 1.31e+00 7.04e+00 +time-adap 28 8.19e-02 1.83e-01 4.12e-01 5.12e+00 +solve 42 1.44e-01 3.74e-01 8.08e-01 1.57e+01 +output 28 8.17e-02 1.82e-01 4.08e-01 5.09e+00 +assemble 42 1.87e-01 3.80e-01 7.86e-01 1.60e+01 +update-ghost 73 1.32e-05 3.49e-04 1.48e-03 2.55e-02 +setup 1 1.99e-01 1.99e-01 1.99e-01 1.99e-01 +--------------------------------------------------------------- +Max number of DoFs: 82891 + +~ RUN 4 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 12 2.59e-01 4.94e-01 1.09e+00 5.93e+00 +time-adap 28 9.71e-02 1.69e-01 3.34e-01 4.72e+00 +solve 42 1.83e-01 3.21e-01 5.57e-01 1.35e+01 +output 28 9.64e-02 1.68e-01 3.32e-01 4.70e+00 +assemble 42 2.02e-01 3.47e-01 8.66e-01 1.46e+01 +update-ghost 73 9.91e-06 3.21e-04 2.67e-03 2.34e-02 +setup 1 2.53e-01 2.53e-01 2.53e-01 2.53e-01 +--------------------------------------------------------------- +Max number of DoFs: 82891 + +~ RUN 5 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 12 2.54e-01 4.71e-01 9.76e-01 5.66e+00 +time-adap 28 8.62e-02 1.51e-01 2.45e-01 4.24e+00 +solve 42 1.31e-01 2.93e-01 5.24e-01 1.23e+01 +output 28 8.58e-02 1.51e-01 2.44e-01 4.21e+00 +assemble 42 1.85e-01 3.21e-01 5.62e-01 1.35e+01 +update-ghost 73 8.89e-06 2.99e-04 1.22e-03 2.18e-02 +setup 1 2.29e-01 2.29e-01 2.29e-01 2.29e-01 +--------------------------------------------------------------- +Max number of DoFs: 82891 + + +======================================================================= + Number of processes = 2 +======================================================================= + +~ RUN 1 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 24 1.63e-01 1.47e-01 6.45e-01 3.53e+00 +time-adap 56 5.22e-02 4.68e-02 1.45e-01 2.62e+00 +solve 84 1.44e-01 1.24e-01 3.70e-01 1.04e+01 +output 56 5.10e-02 4.58e-02 1.42e-01 2.57e+00 +assemble 84 9.23e-02 8.35e-02 2.88e-01 7.02e+00 +update-ghost 146 6.02e-06 6.20e-04 1.66e-02 9.06e-02 +setup 2 1.67e-01 8.37e-02 1.67e-01 1.67e-01 +--------------------------------------------------------------- +Max number of DoFs: 74291 + +~ RUN 2 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 24 1.54e-01 1.67e-01 9.38e-01 4.02e+00 +time-adap 56 5.24e-02 4.53e-02 1.68e-01 2.54e+00 +solve 84 1.21e-01 1.26e-01 3.93e-01 1.06e+01 +output 56 5.19e-02 4.45e-02 1.67e-01 2.49e+00 +assemble 84 9.72e-02 9.17e-02 4.04e-01 7.70e+00 +update-ghost 146 6.87e-06 9.72e-04 3.02e-02 1.42e-01 +setup 2 1.54e-01 7.73e-02 1.55e-01 1.55e-01 +--------------------------------------------------------------- +Max number of DoFs: 74291 + +~ RUN 3 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 24 1.61e-01 1.50e-01 7.33e-01 3.60e+00 +time-adap 56 5.13e-02 4.92e-02 1.72e-01 2.75e+00 +solve 84 1.46e-01 1.41e-01 6.89e-01 1.19e+01 +output 56 5.06e-02 4.82e-02 1.70e-01 2.70e+00 +assemble 84 9.76e-02 9.21e-02 4.44e-01 7.73e+00 +update-ghost 146 4.44e-06 5.52e-04 1.43e-02 8.05e-02 +setup 2 1.50e-01 7.52e-02 1.50e-01 1.50e-01 +--------------------------------------------------------------- +Max number of DoFs: 74291 + +~ RUN 4 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 24 1.50e-01 1.60e-01 7.26e-01 3.83e+00 +time-adap 56 5.64e-02 4.99e-02 1.74e-01 2.80e+00 +solve 84 1.56e-01 1.33e-01 5.21e-01 1.12e+01 +output 56 5.46e-02 4.90e-02 1.71e-01 2.74e+00 +assemble 84 9.72e-02 8.94e-02 2.82e-01 7.51e+00 +update-ghost 146 5.77e-06 7.48e-04 2.03e-02 1.09e-01 +setup 2 1.44e-01 7.22e-02 1.44e-01 1.44e-01 +--------------------------------------------------------------- +Max number of DoFs: 74291 + +~ RUN 5 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 24 1.66e-01 1.53e-01 6.56e-01 3.66e+00 +time-adap 56 5.00e-02 4.68e-02 1.87e-01 2.62e+00 +solve 84 1.61e-01 1.31e-01 4.34e-01 1.10e+01 +output 56 4.92e-02 4.62e-02 1.82e-01 2.59e+00 +assemble 84 1.02e-01 8.85e-02 3.18e-01 7.43e+00 +update-ghost 146 6.89e-06 6.58e-04 1.14e-02 9.60e-02 +setup 2 1.71e-01 8.54e-02 1.71e-01 1.71e-01 +--------------------------------------------------------------- +Max number of DoFs: 74291 + +======================================================================= + Number of processes = 3 +======================================================================= + +~ RUN 1 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 39 1.25e-01 8.52e-02 6.88e-01 3.32e+00 +time-adap 84 4.62e-02 2.62e-02 1.20e-01 2.20e+00 +solve 126 1.59e-01 8.63e-02 4.86e-01 1.09e+01 +output 84 4.45e-02 2.58e-02 1.17e-01 2.16e+00 +assemble 126 7.30e-02 4.57e-02 2.73e-01 5.76e+00 +update-ghost 219 3.52e-06 5.89e-04 1.93e-02 1.29e-01 +setup 3 1.56e-01 5.19e-02 1.56e-01 1.56e-01 +--------------------------------------------------------------- +Max number of DoFs: 72228 + +~ RUN 2 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 39 1.06e-01 8.21e-02 6.14e-01 3.20e+00 +time-adap 84 4.74e-02 2.94e-02 1.32e-01 2.47e+00 +solve 126 1.72e-01 9.58e-02 4.62e-01 1.21e+01 +output 84 4.66e-02 2.89e-02 1.24e-01 2.43e+00 +assemble 126 8.18e-02 4.65e-02 2.57e-01 5.86e+00 +update-ghost 219 3.22e-06 9.72e-04 3.49e-02 2.13e-01 +setup 3 1.86e-01 6.21e-02 1.86e-01 1.86e-01 +--------------------------------------------------------------- +Max number of DoFs: 72228 + +~ RUN 3 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 39 1.16e-01 8.12e-02 6.02e-01 3.17e+00 +time-adap 84 4.54e-02 2.67e-02 1.36e-01 2.24e+00 +solve 126 1.41e-01 8.80e-02 4.69e-01 1.11e+01 +output 84 4.43e-02 2.62e-02 1.35e-01 2.20e+00 +assemble 126 7.10e-02 4.76e-02 2.59e-01 6.00e+00 +update-ghost 219 2.80e-06 7.64e-04 1.89e-02 1.67e-01 +setup 3 1.43e-01 4.78e-02 1.43e-01 1.43e-01 +--------------------------------------------------------------- +Max number of DoFs: 72228 + +~ RUN 4 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 39 1.13e-01 8.56e-02 8.66e-01 3.34e+00 +time-adap 84 4.62e-02 2.62e-02 1.19e-01 2.20e+00 +solve 126 1.58e-01 8.47e-02 5.39e-01 1.07e+01 +output 84 4.43e-02 2.57e-02 1.16e-01 2.16e+00 +assemble 126 7.93e-02 4.61e-02 2.65e-01 5.81e+00 +update-ghost 219 3.20e-06 6.52e-04 1.92e-02 1.43e-01 +setup 3 1.62e-01 5.41e-02 1.62e-01 1.62e-01 +--------------------------------------------------------------- +Max number of DoFs: 72228 + +~ RUN 5 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 39 1.19e-01 8.58e-02 7.04e-01 3.34e+00 +time-adap 84 4.62e-02 2.54e-02 1.39e-01 2.14e+00 +solve 126 1.54e-01 8.86e-02 5.57e-01 1.12e+01 +output 84 4.48e-02 2.50e-02 1.38e-01 2.10e+00 +assemble 126 8.17e-02 4.51e-02 2.23e-01 5.68e+00 +update-ghost 219 2.98e-06 7.45e-04 2.22e-02 1.63e-01 +setup 3 1.84e-01 6.13e-02 1.84e-01 1.84e-01 +--------------------------------------------------------------- +Max number of DoFs: 72228 + +======================================================================= + Number of processes = 4 +======================================================================= + +~ RUN 1 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 56 9.87e-02 4.96e-02 5.50e-01 2.78e+00 +time-adap 112 4.68e-02 1.80e-02 1.33e-01 2.01e+00 +solve 168 1.53e-01 5.65e-02 3.44e-01 9.50e+00 +output 112 4.54e-02 1.76e-02 1.32e-01 1.97e+00 +assemble 168 6.93e-02 2.59e-02 1.67e-01 4.35e+00 +update-ghost 292 2.67e-06 5.14e-04 2.75e-02 1.50e-01 +setup 4 1.32e-01 3.30e-02 1.32e-01 1.32e-01 +--------------------------------------------------------------- +Max number of DoFs: 69285 + +~ RUN 2 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 56 1.03e-01 4.92e-02 4.87e-01 2.75e+00 +time-adap 112 3.79e-02 1.80e-02 1.12e-01 2.02e+00 +solve 168 1.34e-01 5.96e-02 4.30e-01 1.00e+01 +output 112 3.68e-02 1.76e-02 1.10e-01 1.98e+00 +assemble 168 5.41e-02 2.77e-02 1.54e-01 4.66e+00 +update-ghost 292 2.62e-06 4.84e-04 2.16e-02 1.41e-01 +setup 4 1.35e-01 3.38e-02 1.35e-01 1.35e-01 +--------------------------------------------------------------- +Max number of DoFs: 69285 + +~ RUN 3 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 56 9.35e-02 4.80e-02 4.93e-01 2.69e+00 +time-adap 112 4.27e-02 1.61e-02 9.24e-02 1.81e+00 +solve 168 1.20e-01 5.59e-02 3.70e-01 9.39e+00 +output 112 4.18e-02 1.58e-02 9.10e-02 1.77e+00 +assemble 168 6.56e-02 2.61e-02 1.71e-01 4.38e+00 +update-ghost 292 2.88e-06 3.97e-04 1.39e-02 1.16e-01 +setup 4 1.58e-01 3.94e-02 1.58e-01 1.58e-01 +--------------------------------------------------------------- +Max number of DoFs: 69285 + +~ RUN 4 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 56 9.57e-02 4.67e-02 5.04e-01 2.61e+00 +time-adap 112 3.95e-02 1.71e-02 1.14e-01 1.91e+00 +solve 168 1.29e-01 5.97e-02 3.96e-01 1.00e+01 +output 112 3.81e-02 1.68e-02 1.09e-01 1.88e+00 +assemble 168 7.09e-02 2.70e-02 1.62e-01 4.54e+00 +update-ghost 292 2.38e-06 4.21e-04 1.32e-02 1.23e-01 +setup 4 1.18e-01 2.94e-02 1.18e-01 1.18e-01 +--------------------------------------------------------------- +Max number of DoFs: 69285 + +~ RUN 5 ~ + +Adaptivity Results +Number of steps: 42 +Elapsed Time: +--------------------------------------------------------------- +Timer Calls Min(s) Avg(s) Max(s) Total(s) +--------------------------------------------------------------- +refine 56 9.49e-02 5.39e-02 4.96e-01 3.02e+00 +time-adap 112 4.00e-02 1.77e-02 1.39e-01 1.99e+00 +solve 168 1.36e-01 5.97e-02 4.13e-01 1.00e+01 +output 112 3.90e-02 1.74e-02 1.35e-01 1.95e+00 +assemble 168 6.39e-02 2.66e-02 1.78e-01 4.47e+00 +update-ghost 292 2.29e-06 4.50e-04 1.27e-02 1.31e-01 +setup 4 1.54e-01 3.86e-02 1.55e-01 1.55e-01 +--------------------------------------------------------------- +Max number of DoFs: 69285 \ No newline at end of file diff --git a/experiments_results/time_convergence.txt b/experiments_results/time_convergence.txt new file mode 100644 index 000000000..c87fb2ebd --- /dev/null +++ b/experiments_results/time_convergence.txt @@ -0,0 +1,19 @@ +0.02 against 0.04 +L2_error = 4.03e-05 + +0.04 against 0.06 +Computing L2 against baseline +L2_error = 4.76e-04 + +0.06 against 0.08 +Computing L2 against baseline +L2_error = 3.81e-04 + +0.08 against 0.1 +Computing L2 against baseline +L2_error = 1.04e-03 + + +0.01 vs 0.02 +Computing L2 against baseline +L2_error = 9.86e-06 diff --git a/report/Space_Time_adaptivity_report.pdf b/report/Space_Time_adaptivity_report.pdf new file mode 100644 index 000000000..9c5f0b711 Binary files /dev/null and b/report/Space_Time_adaptivity_report.pdf differ diff --git a/src/adaptive/STA_Heat.cpp b/src/adaptive/STA_Heat.cpp new file mode 100644 index 000000000..be1481683 --- /dev/null +++ b/src/adaptive/STA_Heat.cpp @@ -0,0 +1,472 @@ +#include "STA_Heat.hpp" + +void AdaptiveHeat::setup() +{ + + // Create the initial mesh + { + + GridGenerator::hyper_cube(mesh); + //refining level + mesh.refine_global(5); + + } + + + // Initialize the finite element space. + { + fe = std::make_unique>(r); + + quadrature = std::make_unique>(r + 1); + + } + + setup_system(); +} + +void AdaptiveHeat::setup_system(){ + // Initialize the DoF handler. + { + + dof_handler.reinit(mesh); + dof_handler.distribute_dofs(*fe); + + max_dofs = std::max(max_dofs, dof_handler.n_dofs()); + } + + // Initialize the linear system. + { + + const IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs(); + const IndexSet locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler); + + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + // deal.II tutorial 40 also uses SparsityTools::distribute_sparsity_pattern() + // but TrilinosWrappers::SparsityPattern already creates a parallel sparsity pattern + TrilinosWrappers::SparsityPattern sparsity(locally_owned_dofs, + MPI_COMM_WORLD); + // keep_constrained_dofs=false -> Do not treat constrained DoFs as independent unknowns. + // Instead, the constraints are substituted into the matrix structure. + DoFTools::make_sparsity_pattern(dof_handler, + sparsity, + constraints, + /*keep_constrained_dofs=*/false); + sparsity.compress(); + + system_matrix.reinit(sparsity); + + system_rhs.reinit(locally_owned_dofs, MPI_COMM_WORLD); + solution_owned.reinit(locally_owned_dofs, MPI_COMM_WORLD); + solution.reinit(locally_owned_dofs, locally_relevant_dofs, MPI_COMM_WORLD); + } +} + +void AdaptiveHeat::assemble() +{ + // Number of local DoFs for each element. + const unsigned int dofs_per_cell = fe->dofs_per_cell; + + // Number of quadrature points for each element. + const unsigned int n_q = quadrature->size(); + + FEValues fe_values(*fe, + *quadrature, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + // Local matrix and vector. + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector dof_indices(dofs_per_cell); + + // Reset the global matrix and vector, just in case. + system_matrix = 0.0; + system_rhs = 0.0; + + // Evaluation of the old solution on quadrature nodes of current cell. + std::vector solution_old_values(n_q); + + // Evaluation of the gradient of the old solution on quadrature nodes of + // current cell. + std::vector> solution_old_grads(n_q); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (!cell->is_locally_owned()) + continue; + + fe_values.reinit(cell); + + cell_matrix = 0.0; + cell_rhs = 0.0; + + // Evaluate the old solution and its gradient on quadrature nodes. + fe_values.get_function_values(solution, solution_old_values); + fe_values.get_function_gradients(solution, solution_old_grads); + + // cached values useful in the loop + const double inv_dt = 1.0 / delta_t; + const double one_minus_theta = 1.0 - theta; + + + for (unsigned int q = 0; q < n_q; ++q) + { + + const double JxW_q = fe_values.JxW(q); + const Point qp = fe_values.quadrature_point(q); + + const double mu_loc = mu(qp); + + const double f_old_loc = + f(qp, time - delta_t); + const double f_new_loc = f(qp, time); + + + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + // Time derivative. + cell_matrix(i, j) += (inv_dt) * // + fe_values.shape_value(i, q) * // + fe_values.shape_value(j, q) * // + JxW_q; + + // Diffusion. + cell_matrix(i, j) += + theta * mu_loc * // + scalar_product(fe_values.shape_grad(i, q), // + fe_values.shape_grad(j, q)) * // + JxW_q; + } + + // Time derivative. + cell_rhs(i) += (inv_dt) * // + fe_values.shape_value(i, q) * // + solution_old_values[q] * // + JxW_q; + + // Diffusion. + cell_rhs(i) -= (one_minus_theta) * mu_loc * // + scalar_product(fe_values.shape_grad(i, q), // + solution_old_grads[q]) * // + JxW_q; + + // Forcing term. + cell_rhs(i) += + (theta * f_new_loc + (one_minus_theta) * f_old_loc) * // + fe_values.shape_value(i, q) * // + JxW_q; + } + } + + cell->get_dof_indices(dof_indices); + + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + dof_indices, + system_matrix, + system_rhs); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + + // Homogeneous Neumann boundary conditions: we do nothing. +} + +void AdaptiveHeat::solve_time_step() +{ + + TrilinosWrappers::PreconditionAMG preconditioner; + preconditioner.initialize(system_matrix); + + //ReductionControl is a more flexible SolverControl extension + ReductionControl solver_control(/* maxiter = */ 10000, + /* tolerance = */ 1.0e-12, + /* reduce = */ 1.0e-6); + + SolverCG solver(solver_control); + + solution_owned.reinit(dof_handler.locally_owned_dofs(), MPI_COMM_WORLD); + + solver.solve(system_matrix, solution_owned, system_rhs, preconditioner); + + constraints.distribute(solution_owned); + +} + +bool AdaptiveHeat::refine_grid(const unsigned int min_grid_level, + const unsigned int max_grid_level) +{ + Vector estimated_error_per_cell(mesh.n_active_cells()); + + KellyErrorEstimator::estimate( + dof_handler, + QGauss(fe->degree + 1), + std::map *>(), + solution, + estimated_error_per_cell); + + double eta_norm = estimated_error_per_cell.l2_norm(); + + double global_eta; + MPI_Allreduce(&eta_norm, + &global_eta, + 1, + MPI_DOUBLE, + MPI_MAX, + MPI_COMM_WORLD); + + + + if(global_eta < spatial_tol){ + /* do not refine the grid */ + return false; + } + + /* changes to make refine and coarsen fraction dynamic */ + const double error_ratio = eta_norm / spatial_tol; + + // how does it work + // If error_ratio = 2, error is 2x the tol → use higher refine_fraction + // if error_ratio = 1.1, error is a little higher than tol → use lower refine_fraction + const double refine_fraction = std::clamp(0.10 * error_ratio, 0.05, 0.30); //refine_fraction = 0.10 * error_ratio, 0.05 is the min and 0.30 is the max + const double coarsen_fraction = std::clamp(0.01 / error_ratio, 0.005, 0.03); + + + parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction( + mesh, + estimated_error_per_cell, + refine_fraction, + coarsen_fraction, + VectorTools::L1_norm); + /** + * EXPLANATION parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction fuction fields + * 0.3 and 0.03 values explaination (I think it can be useful) + * There is an estimated error for each cell: + * - error is high --> the cell needs more detail + * - error is low --> the cell is already fine + * So: + * - 0.3(refine fraction) --> the worst 30% of cells are split into smaller cells + * - 0.03(coarsen fraction) --> the best 3% of cells are merged into bigger cells + * + * This are default values. + ** How to change them in the future(keep coarsen much smaller than refine (often 5-10x smaller)): + * - if error decreases too slowly, increase refine fraction (0.4 max-ish). + * - if mesh oscillates (refine/coarsen flip-flop), reduce coarsen fraction. + * - if DoFs explode, reduce refine fraction and/or add max refinement level. + */ + + + + // Enforce max and min/max refinement levels + if (mesh.n_levels() > max_grid_level){ + for (const auto &cell : mesh.active_cell_iterators_on_level(max_grid_level)){ + cell->clear_refine_flag(); + } + } + + for (const auto &cell : mesh.active_cell_iterators_on_level(min_grid_level)){ + cell->clear_coarsen_flag(); + } + + + // Prepare solution transfer + parallel::distributed::SolutionTransfer solution_transfer(dof_handler); + // previous_solution is the solution to transfer (tutorial-26) + TrilinosWrappers::MPI::Vector previous_solution(solution); + + mesh.prepare_coarsening_and_refinement(); + solution_transfer.prepare_for_coarsening_and_refinement(previous_solution); + + // Execute refinement + mesh.execute_coarsening_and_refinement(); + + // Rebuild DoFs/matrix/vectors on new mesh + setup_system(); + + // Interpolate old solution onto the new DoF space and apply hanging-node and Dirichlet constraints for continuity. + solution_transfer.interpolate(solution_owned); + constraints.distribute(solution_owned); + + // Update ghosted vector + solution = solution_owned; + return true; +} + +void AdaptiveHeat::output() const +{ + DataOut data_out; + + data_out.add_data_vector(dof_handler, solution, "solution"); + + // Add vector for parallel partition. + std::vector partition_int(mesh.n_active_cells()); + GridTools::get_subdomain_association(mesh, partition_int); + const Vector partitioning(partition_int.begin(), partition_int.end()); + data_out.add_data_vector(partitioning, "partitioning"); + + data_out.build_patches(); + + const std::string output_file_name = "output-mesh"; + + data_out.write_vtu_with_pvtu_record(/* folder = */ "./", + /* basename = */ output_file_name, + /* index = */ timestep_number, + MPI_COMM_WORLD); +} + +void AdaptiveHeat::run() +{ + const unsigned int initial_global_refinement = 6; + + //Checkpoint for rollback + TrilinosWrappers::MPI::Vector old_solution; + + // Setup initial conditions. + { + profiler.tic("setup"); + setup(); + profiler.toc("setup"); + + VectorTools::interpolate(dof_handler, Functions::ZeroFunction(), solution_owned); + solution = solution_owned; + + + old_solution.reinit(solution_owned); + + time = 0.00; + timestep_number = 1; + + + + // Output initial condition. + output(); + } + + // for time error + double tol=1e-2; + double dt_min = 1e-4; + double dt_max = 1e-1; + + + // Time-stepping loop. + while (time < T - 0.5 * delta_t) + { + num_of_steps++; + + profiler.tic("update-ghost"); + //Saving previous solution for error estimation + old_solution = solution_owned; + profiler.toc("update-ghost"); + + + double t_old = time; + double t_attempt = time + delta_t; + time = t_attempt; + + profiler.tic("assemble"); + assemble(); + profiler.toc("assemble"); + + profiler.tic("solve"); + solve_time_step(); + profiler.toc("solve"); + + profiler.tic("time-adap"); + // Error estimation + TrilinosWrappers::MPI::Vector diff = solution_owned; + diff.add(-1.0, old_solution); + + double delta_U = diff.linfty_norm(); + + + // danger for by 0 division + double denom = std::max(delta_U, 1e-14); + //Factor clamp + double factor = std::max(0.3, std::min(tol / denom, 2.0)); + + // Adaptive Rollback + if (delta_U > tol) { + //Rejected step + time =t_old; + + delta_t = 0.9 * delta_t * factor; + //Delta_t clamp + delta_t = std::max(dt_min, std::min(delta_t, dt_max)); + + solution_owned = old_solution; // Back to previous solution + + continue; + } + else { + //Accepted step + time = t_attempt; + timestep_number++; + + delta_t = 0.9 * delta_t * factor; + //Delta_t clamp + delta_t = std::max(dt_min, std::min(delta_t, dt_max)); + + profiler.tic("update-ghost"); + solution = solution_owned; + profiler.toc("update-ghost"); + + profiler.tic("output"); + output(); + profiler.toc("output"); + } + profiler.toc("time-adap"); + + if (time < T - 0.5 * delta_t && timestep_number - last_refine_step >= min_steps_between_refine) + { + + profiler.tic("refine"); + if(refine_grid(initial_global_refinement, initial_global_refinement + 2)){ + last_refine_step = timestep_number; // update last refinement + old_solution.reinit(solution_owned); + + profiler.tic("update-ghost"); + old_solution = solution_owned; + profiler.toc("update-ghost"); + //delta_t = 0.5 * delta_t; + } + profiler.toc("refine"); + + } + } + +} + +double AdaptiveHeat::l2_against_base(const Function & baseline_function){ + + Vector error_per_cell_L2(mesh.n_active_cells()); + + const QGauss quadrature_error(fe->degree + 2); + + VectorTools::integrate_difference( + dof_handler, + solution, + baseline_function, + error_per_cell_L2, + quadrature_error, + VectorTools::L2_norm); + + return error_per_cell_L2.l2_norm(); +} + +void AdaptiveHeat::print_results(){ + pcout << "\nAdaptivity Results" << std::endl; + pcout << "Number of steps: " << num_of_steps << std::endl + << "Elapsed Time: " << std::endl; + profiler.report(MPI_COMM_WORLD); + pcout << "Max number of DoFs: " << max_dofs << std::endl; +} diff --git a/src/adaptive/STA_Heat.hpp b/src/adaptive/STA_Heat.hpp new file mode 100644 index 000000000..d26493f9f --- /dev/null +++ b/src/adaptive/STA_Heat.hpp @@ -0,0 +1,187 @@ +#ifndef STA_HEAT +#define STA_HEAT + +#include +#include +#include +#include // + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // + +#include +#include // +#include +#include +#include // +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include // +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "../timing/Timing.hpp" + +using namespace dealii; + +/** + * Class managing the Heat problem. + */ +class AdaptiveHeat +{ +public: + // Physical dimension (1D, 2D, 3D) + static constexpr unsigned int dim = 3; + + // Class Constructor + AdaptiveHeat(const unsigned int &r_, + const double &T_, + const double &theta_, + const double &delta_t_, + const std::function &)> &mu_, // in this project 1.0 + const std::function &, const double &)> &f_) + : r(r_) + , T(T_) + , theta(theta_) + , delta_t(delta_t_) + , mu(mu_) + , f(f_) + , mpi_size(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)) + , mpi_rank(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)) + , mesh(MPI_COMM_WORLD) + , pcout(std::cout, mpi_rank == 0) + {} + + // Run the simulation for the Heat problem + void run(); + + // Computer L2 error against a baseline solution + double l2_against_base(const Function &); + + void print_results(); + +protected: + // Initialize the mesh and the FE space + void setup(); + + // Initialize the DoFs handler and the system + void setup_system(); + + // Assemble the system + void assemble(); + + // Solve the time step linear system + void solve_time_step(); + + // Apply adaptive mesh refinement + bool refine_grid(const unsigned int, const unsigned int); + + // Output + void output() const; + + // Polynomial degree + const unsigned int r; + + // Final time + const double T; + + // Theta parameter for theta method + const double theta; + + // Timestep + double delta_t; + + // Current time + double time = 0.0; + + // Current timestep number + unsigned int timestep_number = 0; + + // For spatial refinement + double spatial_tol = 5.0e-4; // start point, we can refine + unsigned int min_steps_between_refine = 5; // to not refine too often + unsigned int last_refine_step = 0; // when was the last refinement + + // Diffusion coefficient + std::function &)> mu; + + // Forcing term + std::function &, const double &)> f; + + // Number of MPI processes + const unsigned int mpi_size; + + // Rank of the current MPI process + const unsigned int mpi_rank; + + // Number of steps + unsigned int num_of_steps = 0; + unsigned int max_dofs = 0; + + // Collect data about execution + Profiler profiler; + + // Note: deal.II tutorial uses parallel::distributed::Triangulation. + // fullydistributed is newer and more scalable (does not replicate coarse mesh + // in all mpi processes) however it is more difficult to implement Adaptive Mesh Refinement + // since parallel::distributed::Triangulation already supports many useful features. + // Triangulation + parallel::distributed::Triangulation mesh; + + // Finite element space + std::unique_ptr> fe; + + // Quadrature formula + std::unique_ptr> quadrature; + + // DoF handler + DoFHandler dof_handler; + + // Holds a list of constraints to hold the hanging nodes and the boundary conditions. + AffineConstraints constraints; + + // System matrix + TrilinosWrappers::SparseMatrix system_matrix; + + // System right-hand side + TrilinosWrappers::MPI::Vector system_rhs; + + // System solution, no ghost elements + TrilinosWrappers::MPI::Vector solution_owned; + + // System solution, with ghost elements + TrilinosWrappers::MPI::Vector solution; + + // Output stream, only for process 0 + ConditionalOStream pcout; +}; + +#endif \ No newline at end of file diff --git a/src/adaptive/main.cpp b/src/adaptive/main.cpp new file mode 100644 index 000000000..ce5e50e0d --- /dev/null +++ b/src/adaptive/main.cpp @@ -0,0 +1,95 @@ +#include "STA_Heat.hpp" +//#define COMPARE_WITH_BASE + +#ifdef COMPARE_WITH_BASE +#include "../homogeneous/H_Heat.hpp" +#include +#include +#endif + +int +main(int argc, char *argv[]) +{ + constexpr unsigned int dim = AdaptiveHeat::dim; + + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv); + + const auto mu = [](const Point & /*p*/) { return 1.0; }; + double a = 1.5; + int N = 3; + Point x0(0, 0, 0); + double sigma = 0.5; + const auto g = [&]( const double &t) { + return (std::exp(-a * (std::cos(2*N*M_PI*t) + 1))); + }; + const auto h = [&]( const Point &p) { + return std::exp(-((p-x0)*(p-x0)/std::pow(sigma,2))); + }; + + const auto f = [&](const Point &p, const double &t) { + return g(t)*h(p); + }; + + +#ifdef COMPARE_WITH_BASE + + if(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) > 1){ + if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + std::cout << "Baseline comparing only runs with one process for simplicity." << std::endl + << "Either run with one process or disable comparing with baseline" << std::endl; + return 0; + } + std::unique_ptr> baseline_function; + Heat baseline_heat(/*output_filename = */ "homogeneous.msh", + /* degree = */ 1, + /* T = */ 1.0, + /* theta = */ 0.5, + /* delta_t = */ 0.025, // use a small timestep + mu, + f); + + Vector baseline_serial_solution; + + baseline_heat.run(); + + baseline_serial_solution.reinit(baseline_heat.get_dof_handler().n_dofs()); + baseline_serial_solution = baseline_heat.get_serial_solution(); + + MappingQ1 mapping; + baseline_function = std::make_unique> ( + baseline_heat.get_dof_handler(), + baseline_serial_solution, + mapping + ); + baseline_function->set_time(1.0); + +#endif + + AdaptiveHeat problem(/* degree = */ 1, + /* T = */ 1.0, + /* theta = */ 0.5, + /* delta_t = */ 0.05, + mu, + f); + + problem.run(); + + +#ifdef COMPARE_WITH_BASE + + if(baseline_function){ + baseline_heat.print_results(); + problem.print_results(); + + double L2_err = 0; + + if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0){ + L2_err = problem.l2_against_base(*baseline_function); + } + MPI_Bcast(&L2_err, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + } + +#endif + + return 0; +} \ No newline at end of file diff --git a/src/homogeneous/H_Heat.cpp b/src/homogeneous/H_Heat.cpp new file mode 100644 index 000000000..f3c0196d8 --- /dev/null +++ b/src/homogeneous/H_Heat.cpp @@ -0,0 +1,269 @@ +#include "H_Heat.hpp" + +void +Heat::setup() +{ + + // Create the mesh. + { + + GridGenerator::hyper_cube(mesh); + //refining level + mesh.refine_global(5); + + } + + + // Initialize the finite element space. + { + + fe = std::make_unique>(r); + + quadrature = std::make_unique>(r + 1); + + } + + + // Initialize the DoF handler. + { + + dof_handler.reinit(mesh); + dof_handler.distribute_dofs(*fe); + + } + + // Initialize the linear system. + { + const IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs(); + const IndexSet locally_relevant_dofs = + DoFTools::extract_locally_relevant_dofs(dof_handler); + + TrilinosWrappers::SparsityPattern sparsity(locally_owned_dofs, + MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern(dof_handler, sparsity); + sparsity.compress(); + + system_matrix.reinit(sparsity); + + system_rhs.reinit(locally_owned_dofs, MPI_COMM_WORLD); + solution_owned.reinit(locally_owned_dofs, MPI_COMM_WORLD); + solution.reinit(locally_owned_dofs, locally_relevant_dofs, MPI_COMM_WORLD); + } +} + +void +Heat::assemble() +{ + // Number of local DoFs for each element. + const unsigned int dofs_per_cell = fe->dofs_per_cell; + + // Number of quadrature points for each element. + const unsigned int n_q = quadrature->size(); + + FEValues fe_values(*fe, + *quadrature, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + // Local matrix and vector. + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector dof_indices(dofs_per_cell); + + // Reset the global matrix and vector, just in case. + system_matrix = 0.0; + system_rhs = 0.0; + + // Evaluation of the old solution on quadrature nodes of current cell. + std::vector solution_old_values(n_q); + + // Evaluation of the gradient of the old solution on quadrature nodes of + // current cell. + std::vector> solution_old_grads(n_q); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (!cell->is_locally_owned()) + continue; + + fe_values.reinit(cell); + + cell_matrix = 0.0; + cell_rhs = 0.0; + + // Evaluate the old solution and its gradient on quadrature nodes. + fe_values.get_function_values(solution, solution_old_values); + fe_values.get_function_gradients(solution, solution_old_grads); + + for (unsigned int q = 0; q < n_q; ++q) + { + const double mu_loc = mu(fe_values.quadrature_point(q)); + + const double f_old_loc = + f(fe_values.quadrature_point(q), time - delta_t); + const double f_new_loc = f(fe_values.quadrature_point(q), time); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + // Time derivative. + cell_matrix(i, j) += (1.0 / delta_t) * // + fe_values.shape_value(i, q) * // + fe_values.shape_value(j, q) * // + fe_values.JxW(q); + + // Diffusion. + cell_matrix(i, j) += + theta * mu_loc * // + scalar_product(fe_values.shape_grad(i, q), // + fe_values.shape_grad(j, q)) * // + fe_values.JxW(q); + } + + // Time derivative. + cell_rhs(i) += (1.0 / delta_t) * // + fe_values.shape_value(i, q) * // + solution_old_values[q] * // + fe_values.JxW(q); + + // Diffusion. + cell_rhs(i) -= (1.0 - theta) * mu_loc * // + scalar_product(fe_values.shape_grad(i, q), // + solution_old_grads[q]) * // + fe_values.JxW(q); + + // Forcing term. + cell_rhs(i) += + (theta * f_new_loc + (1.0 - theta) * f_old_loc) * // + fe_values.shape_value(i, q) * // + fe_values.JxW(q); + } + } + + cell->get_dof_indices(dof_indices); + + system_matrix.add(dof_indices, cell_matrix); + system_rhs.add(dof_indices, cell_rhs); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + + // Homogeneous Neumann boundary conditions: we do nothing. +} + +void +Heat::solve_linear_system() +{ + TrilinosWrappers::PreconditionSSOR preconditioner; + preconditioner.initialize( + system_matrix, TrilinosWrappers::PreconditionSSOR::AdditionalData(1.0)); + + ReductionControl solver_control(/* maxiter = */ 10000, + /* tolerance = */ 1.0e-16, + /* reduce = */ 1.0e-6); + + SolverCG solver(solver_control); + + solver.solve(system_matrix, solution_owned, system_rhs, preconditioner); +} + +void +Heat::output() const +{ + DataOut data_out; + + data_out.add_data_vector(dof_handler, solution, "solution"); + + // Add vector for parallel partition. + std::vector partition_int(mesh.n_active_cells()); + GridTools::get_subdomain_association(mesh, partition_int); + const Vector partitioning(partition_int.begin(), partition_int.end()); + data_out.add_data_vector(partitioning, "partitioning"); + + data_out.build_patches(); + + const std::filesystem::path mesh_path(output_file_name); + const std::string output_file_name = "output-" + mesh_path.stem().string(); + + data_out.write_vtu_with_pvtu_record(/* folder = */ "./", + /* basename = */ output_file_name, + /* index = */ timestep_number, + MPI_COMM_WORLD); +} + +void +Heat::run() +{ + // Setup initial conditions. + { + profiler.tic("setup"); + setup(); + profiler.toc("setup"); + + VectorTools::interpolate(dof_handler, FunctionU0(), solution_owned); + solution = solution_owned; + + time = 0.0; + timestep_number = 0; + + // Output initial condition. + output(); + } + + // Time-stepping loop. + while (time < T - 0.5 * delta_t) + { + num_of_steps++; + + time += delta_t; + ++timestep_number; + + + profiler.tic("assemble"); + assemble(); + profiler.toc("assemble"); + + profiler.tic("solve"); + solve_linear_system(); + profiler.toc("solve"); + + // Perform parallel communication to update the ghost values of the + // solution vector. + profiler.tic("update-ghost"); + solution = solution_owned; + profiler.toc("update-ghost"); + + profiler.tic("output"); + output(); + profiler.toc("output"); + } +} + +const DoFHandler & +Heat::get_dof_handler() const +{ + return dof_handler; +} + +Vector +Heat::get_serial_solution() const +{ + Vector serial; + serial.reinit(solution.size()); + for(unsigned int i = 0; i < solution.size(); ++i) + { + serial(i) = solution(i); + } + return serial; +} + +void Heat::print_results(){ + pcout << "Baseline Results" << std::endl; + pcout << "Number of steps: " << num_of_steps << std::endl + << "Elapsed Time: " << std::endl; + profiler.report(MPI_COMM_WORLD); + pcout << "Number of DoFs: " << dof_handler.n_dofs() << std::endl; +} \ No newline at end of file diff --git a/src/homogeneous/H_Heat.hpp b/src/homogeneous/H_Heat.hpp new file mode 100644 index 000000000..33e6cda7a --- /dev/null +++ b/src/homogeneous/H_Heat.hpp @@ -0,0 +1,190 @@ +#ifndef HEAT_HPP +#define HEAT_HPP + +#define _USE_MATH_DEFINES +#include +#include + +#include +#include + +#include +#include // +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + + +#include "../timing/Timing.hpp" + +using namespace dealii; + +/** + * Class managing the differential problem. + */ +class Heat +{ +public: + // Physical dimension (1D, 2D, 3D) + static constexpr unsigned int dim = 3; + + // Initial condition. + class FunctionU0 : public Function + { + public: + // Constructor. + FunctionU0() = default; + + // Evaluation of the function. + virtual double + value(const Point &/*p*/, + const unsigned int /*component*/ = 0) const override + { + return 0.0; + } + }; + + // Constructor. + Heat(const std::string &output_file_name_, + const unsigned int &r_, + const double &T_, + const double &theta_, + const double &delta_t_, + const std::function &)> &mu_, + const std::function &, const double &)> &f_) + : output_file_name(output_file_name_) + , r(r_) + , T(T_) + , theta(theta_) + , delta_t(delta_t_) + , mu(mu_) + , f(f_) + , mpi_size(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)) + , mpi_rank(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)) + , mesh(MPI_COMM_WORLD) + , pcout(std::cout, mpi_rank == 0) + {} + + // Run the time-dependent simulation. + void + run(); + + // function needed for comparison with adaptive problem + Vector get_serial_solution() const; + + const DoFHandler& get_dof_handler() const; + + void print_results(); + +protected: + // Initialization. + void + setup(); + + // System assembly. + void + assemble(); + + // System solution. + void + solve_linear_system(); + + // Output. + void + output() const; + + // Name of the mesh. + const std::string output_file_name; + + // Polynomial degree. + const unsigned int r; + + // Final time. + const double T; + + // Theta parameter for the theta method. + const double theta; + + // Time step. + const double delta_t; + + // Current time. + double time = 0.0; + + // Current timestep number. + unsigned int timestep_number = 0; + + // Diffusion coefficient. + std::function &)> mu; + + // Forcing term. + std::function &, const double &)> f; + + // Number of MPI processes. + const unsigned int mpi_size; + + // Rank of the current MPI process. + const unsigned int mpi_rank; + + // Number of steps + unsigned int num_of_steps = 0; + + // Collect data about execution + Profiler profiler; + + // Triangulation. + parallel::distributed::Triangulation mesh; + + // Finite element space. + std::unique_ptr> fe; + + // Quadrature formula. + std::unique_ptr> quadrature; + + // DoF handler. + DoFHandler dof_handler; + + // System matrix. + TrilinosWrappers::SparseMatrix system_matrix; + + // System right-hand side. + TrilinosWrappers::MPI::Vector system_rhs; + + // System solution, without ghost elements. + TrilinosWrappers::MPI::Vector solution_owned; + + // System solution, with ghost elements. + TrilinosWrappers::MPI::Vector solution; + + // Output stream for process 0. + ConditionalOStream pcout; +}; + +#endif \ No newline at end of file diff --git a/src/homogeneous/main.cpp b/src/homogeneous/main.cpp new file mode 100644 index 000000000..066299d64 --- /dev/null +++ b/src/homogeneous/main.cpp @@ -0,0 +1,38 @@ +#include "H_Heat.hpp" + +// Main function. +int +main(int argc, char *argv[]) +{ + constexpr unsigned int dim = Heat::dim; + + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv); + + const auto mu = [](const Point & /*p*/) { return 1.0; }; + double a = 1.5; + int N = 3; + Point x0(0, 0, 0); + double sigma = 0.5; + const auto g = [&]( const double &t) { + return (std::exp(-a * (std::cos(2*N*M_PI*t) + 1))); + }; + const auto h = [&]( const Point &p) { + return std::exp(-((p-x0)*(p-x0)/std::pow(sigma,2))); + }; + + const auto f = [&](const Point &p, const double &t) { + return g(t)*h(p); + }; + + Heat problem(/*ouput_filename = */ "homogeneus.msh", + /* degree = */ 1, + /* T = */ 1.0, + /* theta = */ 0.5, + /* delta_t = */ 0.05, + mu, + f); + + problem.run(); + + return 0; +} \ No newline at end of file diff --git a/src/timing/Timing.cpp b/src/timing/Timing.cpp new file mode 100644 index 000000000..9a02039ec --- /dev/null +++ b/src/timing/Timing.cpp @@ -0,0 +1,79 @@ +#include "Timing.hpp" + +void Profiler::tic(const std::string& name) +{ + timers[name].tic(); +} + +void Profiler::toc(const std::string& name) +{ + timers[name].toc(); +} + +void Profiler::report(MPI_Comm comm) const +{ + int rank; + MPI_Comm_rank(comm, &rank); + + // Only rank 0 prints final table + struct AggTimer { + double total = 0.0; + double min = 0.0; + double max = 0.0; + int calls = 0; + }; + + if (rank == 0) + { + std::cout << "---------------------------------------------------------------\n"; + std::cout << std::left + << std::setw(18) << "Timer" + << std::setw(10) << "Calls" + << std::setw(10) << "Min(s)" + << std::setw(10) << "Avg(s)" + << std::setw(10) << "Max(s)" + << std::setw(10) << "Total(s)" + << "\n"; + std::cout << "---------------------------------------------------------------\n"; + } + + for (const auto& [name, local] : timers) + { + AggTimer global; + + // Reduce TOTAL + MPI_Reduce(&local.total, &global.total, 1, + MPI_DOUBLE, MPI_MAX, 0, comm); + + // Reduce CALLS + MPI_Reduce(&local.calls, &global.calls, 1, + MPI_INT, MPI_SUM, 0, comm); + + // Reduce MIN (best case across ranks) + MPI_Reduce(&local.min, &global.min, 1, + MPI_DOUBLE, MPI_MIN, 0, comm); + + // Reduce MAX (worst case across ranks) + MPI_Reduce(&local.max, &global.max, 1, + MPI_DOUBLE, MPI_MAX, 0, comm); + + if (rank == 0) + { + double avg = global.total / global.calls; + + std::cout << std::left + << std::setw(18) << name + << std::setw(10) << global.calls + << std::setw(10) << global.min + << std::setw(10) << avg + << std::setw(10) << global.max + << std::setw(10) << global.total + << "\n"; + } + } + + if (rank == 0) + { + std::cout << "---------------------------------------------------------------\n"; + } +} \ No newline at end of file diff --git a/src/timing/Timing.hpp b/src/timing/Timing.hpp new file mode 100644 index 000000000..ed1861c0f --- /dev/null +++ b/src/timing/Timing.hpp @@ -0,0 +1,41 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + + +struct Timer { + double total = 0.0; + double start = 0.0; + + int calls = 0; + + double min = std::numeric_limits::max(); + double max = 0.0; + + void tic() { + start = MPI_Wtime(); // Start timer + } + + void toc() { + double end = MPI_Wtime(); + total += end - start; // Update total time + ++calls; + min = std::min(min, end - start); // update min + max = std::max(max, end - start); // update max + } +}; + +class Profiler { +public: + void tic(const std::string& name); + void toc(const std::string& name); + void report(MPI_Comm) const; + +private: + std::unordered_map timers; +}; \ No newline at end of file