Merge branch 'hannes-philipp' of git.gfz-potsdam.de:naaice/tug into hannes-philipp

.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:

.readthedocs.yaml

@@ -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

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)
+add_subdirectory(examples)

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
 

doc/ValidationHetDiff.org

@@ -1,12 +1,12 @@
-#+TITLE: 2D Validation Examples
+#+TITLE: Validation Examples for 2D Heterogeneous Diffusion
 #+AUTHOR: MDL <delucia@gfz-potsdam.de>
-#+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(ReacTran)
-library(deSolve)
+  library(deSolve)
 
-## harmonic mean
+  ## harmonic mean
-harm <- function(x,y) {
+  harm <- function(x,y) {
-  if (length(x) != 1 || length(y) != 1)
+    if (length(x) != 1 || length(y) != 1)
-  stop("x & y have different lengths")
+    stop("x & y have different lengths")
-  2/(1/x+1/y)
+    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])
   }
-  ret
-}
 
-## Construction of the 2D grid
+  N     <- 11          # number of grid cells
-x.grid <- setup.grid.1D(x.up = 0, L = L, N = N)
+  ini   <- 1           # initial value at x=0
-y.grid <- setup.grid.1D(x.up = 0, L = L, N = N)
+  N2    <- ceiling(N/2)
-grid2D <- setup.grid.2D(x.grid, y.grid)
+  L     <- 10          # domain side
-dx <- dy <- L/N
 
-D.grid <- list()
+  ## Define diff.coeff per cell, in  4 quadrants
-## Diffusion coefs on x-interfaces
+  alphas <- matrix(0, N, N) 
-D.grid$x.int <- apply(alphas, 1, cmpharm)
+  alphas[1:N2, 1:N2] <- 1
-## Diffusion coefs on y-interfaces
+  alphas[1:N2, seq(N2+1,N)] <- 0.1
-D.grid$y.int <- t(apply(alphas, 2, cmpharm))
+  alphas[seq(N2+1,N), 1:N2] <- 0.01
+  alphas[seq(N2+1,N), seq(N2+1,N)] <- 0.001
 
-# The model
+  cmpharm <- function(x) {
-Diff2Dc <- function(t, y, parms) {
+    y <- c(0, x, 0)
-  CONC  <- matrix(nrow = N, ncol = N, data = y)
+    ret <- numeric(length(x)+1)
-  dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC
+    for (i in seq(2, length(y))) {
-  return(list(dCONC))
+      ret[i-1] <- harm(y[i], y[i-1])
-}
+    }
+    ret
+  }
 
-## initial condition: 0 everywhere, except in central point
+  ## Construction of the 2D grid
-y <- matrix(nrow = N, ncol = N, data = 0)
+  x.grid <- setup.grid.1D(x.up = 0, L = L, N = N)
-y[N2, N2] <- ini  # initial concentration in the central point...
+  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
+  D.grid <- list()
-times <- 0:10
+  ## Diffusion coefs on x-interfaces
-outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = NULL,
+  D.grid$x.int <- apply(alphas, 1, cmpharm)
-               dim = c(N, N), lrw = 1860000)
+  ## 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
 

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 
+Welcome to the documentation of the TUG project, a simulation program
-generally, for solving the following differential equation
+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

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)

src/BTCSv2.cpp

@@ -278,10 +278,10 @@ static VectorXd EigenLUAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
 static VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
     uint32_t n = b.size();
 
-    VectorXd a_diag(n);
+    Eigen::VectorXd a_diag(n);
-    VectorXd b_diag(n);
+    Eigen::VectorXd b_diag(n);
-    VectorXd c_diag(n);
+    Eigen::VectorXd c_diag(n);
-    VectorXd x_vec = b;
+    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<double> A;
     VectorXd b;
 
-    // const MatrixXd &alphaX = grid.getAlphaX(); // TODO check if this helps performance
     MatrixXd alphaX = grid.getAlphaX();
     MatrixXd alphaY = grid.getAlphaY();
     vector<BoundaryElement> 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<SparseMatrix<double>> solver; // TODO what is this?
+        SparseLU<SparseMatrix<double>> 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) {
+    if (grid.getDim() == 1) {
-        BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads);
-    } else 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) {
+    if (grid.getDim() == 1) {
-        BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads);
-    } else 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!");
     }