diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index e6af65d..c2608e2 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -4,6 +4,7 @@ stages: - build - test - static_analyze + - doc build_release: stage: build @@ -22,6 +23,23 @@ test: script: - ./build/test/testTug +pages: + stage: doc + image: python:slim + before_script: + - apt-get update && apt-get install --no-install-recommends -y graphviz imagemagick doxygen make + - pip install --upgrade pip && pip install Sphinx Pillow breathe sphinx-rtd-theme + - mkdir public + script: + - pushd docs_sphinx + - make html + - popd && mv docs_sphinx/_build/html/* public/ + artifacts: + paths: + - public + rules: + - if: $CI_COMMIT_REF_NAME == $CI_DEFAULT_BRANCH + lint: stage: static_analyze before_script: diff --git a/.readthedocs.yaml b/.readthedocs.yaml deleted file mode 100644 index 3d4131c..0000000 --- a/.readthedocs.yaml +++ /dev/null @@ -1,32 +0,0 @@ -# .readthedocs.yaml -# Read the Docs configuration file -# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details - -# Required -version: 2 - -# Set the OS, Python version and other tools you might need -build: - os: ubuntu-22.04 - tools: - python: "3.11" - # You can also specify other tool versions: - # nodejs: "19" - # rust: "1.64" - # golang: "1.19" - -# Build documentation in the "docs/" directory with Sphinx -sphinx: - configuration: docs_sphinx/conf.py - -# Optionally build your docs in additional formats such as PDF and ePub -formats: - - pdf -# - epub - -# Optional but recommended, declare the Python requirements required -# to build your documentation -# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html -python: - install: - - requirements: docs_sphinx/requirements.txt diff --git a/CMakeLists.txt b/CMakeLists.txt index 9cd2625..a76e995 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,8 @@ set(CMAKE_CXX_STANDARD 17) find_package(Eigen3 REQUIRED NO_MODULE) find_package(OpenMP) -find_package(easy_profiler) +# find_package(easy_profiler) +# option(EASY_OPTION_LOG "Verbose easy_profiler" 1) ## SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -mfma") option(TUG_USE_OPENMP "Compile with OpenMP support" ON) @@ -40,4 +41,4 @@ if(TUG_ENABLE_TESTING) add_subdirectory(test) endif() -add_subdirectory(examples) \ No newline at end of file +add_subdirectory(examples) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 9408089..575249f 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -1,6 +1,7 @@ #+TITLE: Finite Difference Schemes for the numerical solution of heterogeneous diffusion equation in 2D #+LaTeX_CLASS_OPTIONS: [a4paper,10pt] #+LATEX_HEADER: \usepackage{fullpage} +#+LATEX_HEADER: \usepackage{charter} #+LATEX_HEADER: \usepackage{amsmath, systeme, cancel, xcolor} #+OPTIONS: toc:nil diff --git a/doc/ADI_scheme.pdf b/doc/ADI_scheme.pdf index 62f2f9f..1e389d7 100644 Binary files a/doc/ADI_scheme.pdf and b/doc/ADI_scheme.pdf differ diff --git a/doc/ValidationHetDiff.org b/doc/ValidationHetDiff.org index c9b5613..bb46e3a 100644 --- a/doc/ValidationHetDiff.org +++ b/doc/ValidationHetDiff.org @@ -1,12 +1,12 @@ -#+TITLE: 2D Validation Examples +#+TITLE: Validation Examples for 2D Heterogeneous Diffusion #+AUTHOR: MDL -#+DATE: 2023-07-31 +#+DATE: 2023-08-26 #+STARTUP: inlineimages #+LATEX_CLASS_OPTIONS: [a4paper,9pt] #+LATEX_HEADER: \usepackage{fullpage} #+LATEX_HEADER: \usepackage{amsmath, systeme} #+LATEX_HEADER: \usepackage{graphicx} -#+LATEX_HEADER: \usepackage{} +#+LATEX_HEADER: \usepackage{charter} #+OPTIONS: toc:nil @@ -38,7 +38,7 @@ constant in 4 quadrants: The relevant part of the R script used to produce these results is presented in listing 1; the whole script is at [[file:scripts/HetDiff.R]]. A visualization of the output of the reference simulation is given in -figure [[#fig:1][1]]. +figure [[fig:1][1]]. Note: all results from this script are stored in the =outc= matrix by the =deSolve= function. I stored a different version into @@ -47,68 +47,69 @@ for each time step including initial conditions) and 121 rows, one for each domain element, with no headers. #+caption: Result of ReacTran/deSolve solution of the above problem at 4 +#+name: fig:1 [[./images/deSolve_AlphaHet1.png]] #+name: lst:1 #+begin_src R :language R :frame single :caption Listing 1, generate reference simulation using R packages deSolve/ReacTran :captionpos b :label lst:1 -library(ReacTran) -library(deSolve) + library(ReacTran) + library(deSolve) -## harmonic mean -harm <- function(x,y) { - if (length(x) != 1 || length(y) != 1) - stop("x & y have different lengths") - 2/(1/x+1/y) -} - -N <- 11 # number of grid cells -ini <- 1 # initial value at x=0 -N2 <- ceiling(N/2) -L <- 10 # domain side - -## Define diff.coeff per cell, in 4 quadrants -alphas <- matrix(0, N, N) -alphas[1:N2, 1:N2] <- 1 -alphas[1:N2, seq(N2+1,N)] <- 0.1 -alphas[seq(N2+1,N), 1:N2] <- 0.01 -alphas[seq(N2+1,N), seq(N2+1,N)] <- 0.001 - -cmpharm <- function(x) { - y <- c(0, x, 0) - ret <- numeric(length(x)+1) - for (i in seq(2, length(y))) { - ret[i-1] <- harm(y[i], y[i-1]) + ## harmonic mean + harm <- function(x,y) { + if (length(x) != 1 || length(y) != 1) + stop("x & y have different lengths") + 2/(1/x+1/y) } - ret -} -## Construction of the 2D grid -x.grid <- setup.grid.1D(x.up = 0, L = L, N = N) -y.grid <- setup.grid.1D(x.up = 0, L = L, N = N) -grid2D <- setup.grid.2D(x.grid, y.grid) -dx <- dy <- L/N + N <- 11 # number of grid cells + ini <- 1 # initial value at x=0 + N2 <- ceiling(N/2) + L <- 10 # domain side -D.grid <- list() -## Diffusion coefs on x-interfaces -D.grid$x.int <- apply(alphas, 1, cmpharm) -## Diffusion coefs on y-interfaces -D.grid$y.int <- t(apply(alphas, 2, cmpharm)) + ## Define diff.coeff per cell, in 4 quadrants + alphas <- matrix(0, N, N) + alphas[1:N2, 1:N2] <- 1 + alphas[1:N2, seq(N2+1,N)] <- 0.1 + alphas[seq(N2+1,N), 1:N2] <- 0.01 + alphas[seq(N2+1,N), seq(N2+1,N)] <- 0.001 -# The model -Diff2Dc <- function(t, y, parms) { - CONC <- matrix(nrow = N, ncol = N, data = y) - dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC - return(list(dCONC)) -} + cmpharm <- function(x) { + y <- c(0, x, 0) + ret <- numeric(length(x)+1) + for (i in seq(2, length(y))) { + ret[i-1] <- harm(y[i], y[i-1]) + } + ret + } -## initial condition: 0 everywhere, except in central point -y <- matrix(nrow = N, ncol = N, data = 0) -y[N2, N2] <- ini # initial concentration in the central point... + ## Construction of the 2D grid + x.grid <- setup.grid.1D(x.up = 0, L = L, N = N) + y.grid <- setup.grid.1D(x.up = 0, L = L, N = N) + grid2D <- setup.grid.2D(x.grid, y.grid) + dx <- dy <- L/N -## solve for 10 time units -times <- 0:10 -outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = NULL, - dim = c(N, N), lrw = 1860000) + D.grid <- list() + ## Diffusion coefs on x-interfaces + D.grid$x.int <- apply(alphas, 1, cmpharm) + ## Diffusion coefs on y-interfaces + D.grid$y.int <- t(apply(alphas, 2, cmpharm)) + + # The model + Diff2Dc <- function(t, y, parms) { + CONC <- matrix(nrow = N, ncol = N, data = y) + dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC + return(list(dCONC)) + } + + ## initial condition: 0 everywhere, except in central point + y <- matrix(nrow = N, ncol = N, data = 0) + y[N2, N2] <- ini # initial concentration in the central point... + + ## solve for 10 time units + times <- 0:10 + outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = NULL, + dim = c(N, N), lrw = 1860000) #+end_src diff --git a/docs_sphinx/index.rst b/docs_sphinx/index.rst index 740fbbc..6c5d9ce 100644 --- a/docs_sphinx/index.rst +++ b/docs_sphinx/index.rst @@ -5,12 +5,26 @@ Welcome to Tug's documentation! =============================== -Welcome to the documentation of the TUG project, a simulation program -for solving one- and two-dimensional diffusion problems with heterogeneous diffusion coefficients, more -generally, for solving the following differential equation + +Welcome to the documentation of the TUG project, a simulation program +for solving transport equations in one- and two-dimensional uniform +grids using cell centered finite differences. + +--------- +Diffusion +--------- + +TUG can solve diffusion problems with heterogeneous and anisotropic +diffusion coefficients. The partial differential equation expressing +diffusion reads: .. math:: - \frac{\partial C}{\partial t} = \alpha_x \frac{\partial^2 C}{\partial x^2} + \alpha_y \frac{\partial^2 C}{\partial y^2}. + \frac{\partial C}{\partial t} = \nabla \cdot \left[ \mathbf{\alpha} \nabla C \right] + +In 2D, the equation reads: + +.. math:: + \frac{\partial C}{\partial t} = \frac{\partial}{\partial x}\left[ \alpha_x \frac{\partial C}{\partial x}\right] + \frac{\partial}{\partial y}\left[ \alpha_y \frac{\partial C}{\partial y}\right] .. toctree:: :maxdepth: 2 diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 91c998a..a4ebaaf 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -4,6 +4,7 @@ add_executable(BTCS_1D_proto_example BTCS_1D_proto_example.cpp) add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp) add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp) add_executable(profiling_openmp profiling_openmp.cpp) + target_link_libraries(FTCS_1D_proto_example tug) target_link_libraries(FTCS_2D_proto_example tug) target_link_libraries(BTCS_1D_proto_example tug) diff --git a/src/BTCSv2.cpp b/src/BTCSv2.cpp index 7b170c7..49e4d1c 100644 --- a/src/BTCSv2.cpp +++ b/src/BTCSv2.cpp @@ -278,10 +278,10 @@ static VectorXd EigenLUAlgorithm(SparseMatrix &A, VectorXd &b) { static VectorXd ThomasAlgorithm(SparseMatrix &A, VectorXd &b) { uint32_t n = b.size(); - VectorXd a_diag(n); - VectorXd b_diag(n); - VectorXd c_diag(n); - VectorXd x_vec = b; + Eigen::VectorXd a_diag(n); + Eigen::VectorXd b_diag(n); + Eigen::VectorXd c_diag(n); + Eigen::VectorXd x_vec = b; // Fill diagonals vectors b_diag[0] = A.coeff(0, 0); @@ -367,7 +367,6 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solver SparseMatrix A; VectorXd b; - // const MatrixXd &alphaX = grid.getAlphaX(); // TODO check if this helps performance MatrixXd alphaX = grid.getAlphaX(); MatrixXd alphaY = grid.getAlphaY(); vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); @@ -385,7 +384,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solver b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, colMax, i, sx, sy); - SparseLU> solver; // TODO what is this? + SparseLU> solver; row_t1 = solverFunc(A, b); @@ -417,10 +416,10 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solver // entry point for EigenLU solver; differentiate between 1D and 2D grid static void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { - if (grid.getDim() == 2) { - BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads); - } else if (grid.getDim() == 1) { + if (grid.getDim() == 1) { BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); + } else if (grid.getDim() == 2) { + BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads); } else { throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); } @@ -428,10 +427,10 @@ static void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { // entry point for Thomas algorithm solver; differentiate 1D and 2D grid static void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) { - if (grid.getDim() == 2) { - BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads); - } else if (grid.getDim() == 1) { + if (grid.getDim() == 1) { BTCS_1D(grid, bc, timestep, ThomasAlgorithm); + } else if (grid.getDim() == 2) { + BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads); } else { throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); }