diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 5b38c99..2e36801 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -3,20 +3,54 @@ image: sobc/gitlab-ci stages: - build - test + - static_analyze + - dynamic_analyze before_script: - apt-get update && apt-get install -y libeigen3-dev build: stage: build + artifacts: + paths: + - build/app/1D + - build/app/2D + - build/app/Comp2D + expire_in: 100s script: - mkdir build && cd build - cmake .. - make -lint: +run_1D: stage: test + dependencies: + - build script: - - mkdir build && cd build + - ./build/app/1D + +run_2D: + stage: test + dependencies: + - build + script: + - ./build/app/2D + +lint: + stage: static_analyze + script: + - mkdir lint && cd lint - cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*,readability-*, modernize-*" .. - make + +memcheck_1D: + stage: dynamic_analyze + script: + - cd build/app + - valgrind --tool=memcheck --leak-check=full --show-leak-kinds=all --track-origins=yes ./1D 1>/dev/null + +memcheck_2D: + stage: dynamic_analyze + script: + - cd build/app + - valgrind --tool=memcheck --leak-check=full --show-leak-kinds=all --track-origins=yes ./2D 1>/dev/null diff --git a/CMakeLists.txt b/CMakeLists.txt index a18d9cc..068023f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,4 +7,5 @@ set(CMAKE_CXX_STANDARD 14) find_package(Eigen3 REQUIRED NO_MODULE) +add_subdirectory(app) add_subdirectory(src) diff --git a/Comp.R b/Comp.R new file mode 100644 index 0000000..482e022 --- /dev/null +++ b/Comp.R @@ -0,0 +1,227 @@ +## Time-stamp: "Last modified 2022-01-25 11:22:30 delucia" +library(ReacTran) +library(deSolve) +options(width=114) + + +Diffusion <- function (t, Y, parms){ + tran <- tran.1D(C = Y, C.up = 1, C.down = 0, D = parms$D, dx = xgrid) + return(list(tran$dC)) +} + + +## sol 100 +## Set initial conditions +N <- 100 +xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = N) +x <- xgrid$x.mid +D.coeff <- 1E-3 +C0 <- 1 ## Initial concentration (mg/L) +X0 <- 0 ## Location of initial concentration (m) +Yini <- c(C0, rep(0,N-1)) +parms1 <- list(D=D.coeff) +times <- seq(from = 0, to = 5, by = 1) +system.time({ + out100 <- ode.1D(y = Yini, times = times, func = Diffusion, + parms = parms1, dimens = N) +}) + + +## sol 1000 +## Set initial conditions +N <- 1000 +xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = N) +x <- xgrid$x.mid +D.coeff <- 1E-3 +C0 <- 1 ## Initial concentration (mg/L) +X0 <- 0 ## Location of initial concentration (m) +Yini <- c(C0, rep(0,N-1)) +parms1 <- list(D=D.coeff) +times <- seq(from = 0, to = 5, by = 1) + +system.time({ + out1000 <- ode.1D(y = Yini, times = times, func = Diffusion, + parms = parms1, dimens = N) +}) + +## sol 5000 +## Set initial conditions +N <- 5000 +xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = N) +x <- xgrid$x.mid +D.coeff <- 1E-3 +C0 <- 1 ## Initial concentration (mg/L) +X0 <- 0 ## Location of initial concentration (m) +Yini <- c(C0, rep(0,N-1)) +parms1 <- list(D=D.coeff) +times <- seq(from = 0, to = 5, by = 1) +system.time({ + out5000 <- + ode.1D(y = Yini, times = times, func = Diffusion, parms = parms1, + dimens = N) +}) + +## ./mdl_main 100 > Sol100.dat +## ./mdl_main 1000 > Sol1000.dat +## ./mdl_main 5000 > Sol5000.dat +a100 <- data.table::fread("build/src/Sol100.dat") +a1000 <- data.table::fread("build/src/Sol1000.dat") +a5000 <- data.table::fread("build/src/Sol5000.dat") + + +## run: phreeqc pqc.in pqc.out /PATH_TO_PHREEQC.DAT +p100 <- data.table::fread("./selected_output_1.sel") + +pqc <- subset(p100, subset = time==5) |> + subset(subset = soln > 0 & soln < 101) |> + subset(select = Mg) + +cairo_pdf("test100.pdf", family="serif") +matplot.1D(out100, type = "l", subset = time == 4, xaxs="i", xlab="grid element", ylab="C", + main="Comparison ode.1D vs btcs, discretization 100", lwd=2) +lines(a100$inp1, lwd=2, col="blue") +lines(a100$inp2, lwd=2, col="red") +lines(a100$inp3, lwd=2, col="green4") +lines(pqc$Mg, lwd=2, col="orange") +legend("topright", c("R deSolve ode.1D reference, ts=1", + "btcs, ts=1", "btcs, ts=0.1","btcs, ts=5","PHREEQC ts=1"), lty="solid", + lwd=2, col=c("black","blue","red", "green4", "orange"), + bty="n") +dev.off() + +cairo_pdf("test1000.pdf", family="serif") +matplot.1D(out1000, type = "l", subset = time == 4, xaxs="i", xlab="grid element", ylab="C", + main="Comparison ode.1D vs btcs, discretization 1000", lwd=2) +lines(a1000$inp1, lwd=2, col="blue") +lines(a1000$inp2, lwd=2, col="red") +lines(a1000$inp3, lwd=2, col="green4") +legend("topright", c("R deSolve ode.1D reference, ts=1", + "btcs, ts=1", "btcs, ts=0.1","btcs, ts=5"), lty="solid", + lwd=2, col=c("black","blue","red", "green4"), + bty="n") +dev.off() + +cairo_pdf("test5000.pdf", family="serif") +matplot.1D(out5000, type = "l", subset = time == 4, xaxs="i", xlab="grid element", ylab="C", + main="Comparison ode.1D vs btcs, discretization 5000", lwd=2) +lines(a5000$inp1, lwd=2, col="blue") +lines(a5000$inp2, lwd=2, col="red") +lines(a5000$inp3, lwd=2, col="green4") +legend("topright", c("R deSolve ode.1D reference, ts=1", + "btcs, ts=1", "btcs, ts=0.1","btcs, ts=5"), lty="solid", + lwd=2, col=c("black","blue","red", "green4"), + bty="n") +dev.off() + + +############## old stuff + +a <- data.table::fread("build/src/Sol1000.dat") + +erfc <- function(x) 2 * pnorm(x * sqrt(2), lower = FALSE) + + +Analytical <- function(x, t, a, v) { + erfc +} + +cairo_pdf("test.pdf", family="serif") +matplot.1D(out, type = "l", subset = time == 4, xaxs="i", xlab="grid element", ylab="C", + main="Comparison ode.1D vs btcs", lwd=2) +lines(a$inp1, lwd=2, col="blue") +lines(a$inp2, lwd=2, col="red") +lines(a$inp3, lwd=2, col="green4") +legend("topright", c("R deSolve ode.1D reference, ts=1", + "btcs, ts=1", "btcs, ts=0.1","btcs, ts=5"), lty="solid", + lwd=2, col=c("black","blue","red", "green4"), + bty="n") +dev.off() + + +## sol 5000 +## Set initial conditions +N <- 5000 +xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = N) +x <- xgrid$x.mid +D.coeff <- 1E-3 +C0 <- 1 ## Initial concentration (mg/L) +X0 <- 0 ## Location of initial concentration (m) +Yini <- c(C0, rep(0,N-1)) + +parms1 <- list(D=D.coeff) + +times <- seq(from = 0, to = 5, by = 1) +system.time( + out5000 <- ode.1D(y = Yini, times = times, func = Diffusion, + parms = parms1, dimens = N)) + +a <- data.table::fread("build/src/Sol5000.dat") + +cairo_pdf("test.pdf", family="serif") +matplot.1D(out5000, type = "l", subset = time == 4, xaxs="i", xlab="grid element", ylab="C", + main="Comparison ode.1D vs btcs", lwd=2) +lines(a$inp1, lwd=2, col="blue") +lines(a$inp2, lwd=2, col="red") +lines(a$inp3, lwd=2, col="green4") +legend("topright", c("R deSolve ode.1D reference, ts=1", + "btcs, ts=1", "btcs, ts=0.1","btcs, ts=5"), lty="solid", + lwd=2, col=c("black","blue","red", "green4"), + bty="n") +dev.off() + + + +### ## +## diffusR <- function(t, C, parameters) { +## deltax <- c (0.5*delx, rep(delx, numboxes-1), 0.5*delx) +## Flux <- -D*diff(c(0, C, 0))/deltax +## dC <- -diff(Flux)/delx + +## list(dC) # the output +## } + +## ## ================== +## ## Model application +## ## ================== + +## ## the model parameters: +## D <- 0.3 # m2/day diffusion rate +## r <- 0.01 # /day net growth rate +## delx <- 1 # m thickness of boxes +## numboxes <- 50 + +## ## distance of boxes on plant, m, 1 m intervals +## Distance <- seq(from = 0.5, by = delx, length.out = numboxes) + +## ## Initial conditions, ind/m2 +## ## aphids present only on two central boxes +## C <- rep(0, times = numboxes) +## C[30:31] <- 1 +## state <- c(C = C) # initialise state variables + +## ## RUNNING the model: +## times <- seq(0, 200, by = 1) # output wanted at these time intervals +## out <- ode.band(state, times, diffusR, parms = 0, +## nspec = 1, names = "C") + +## ## ================ +## ## Plotting output +## ## ================ +## image(out, grid = Distance, method = "filled.contour", +## xlab = "time, days", ylab = "Distance on plant, m", +## main = "Aphid density on a row of plants") + +## matplot.1D(out, grid = Distance, type = "l", +## subset = time %in% seq(0, 200, by = 10)) + +## # add an observed dataset to 1-D plot (make sure to use correct name): +## data <- cbind(dist = c(0,10, 20, 30, 40, 50, 60), +## Aphid = c(0,0.1,0.25,0.5,0.25,0.1,0)) + +## matplot.1D(out, grid = Distance, type = "l", +## subset = time %in% seq(0, 200, by = 10), +## obs = data, obspar = list(pch = 18, cex = 2, col="red")) + +## ## Not run: +## plot.1D(out, grid = Distance, type = "l") + diff --git a/Comp2D.R b/Comp2D.R new file mode 100644 index 0000000..7111103 --- /dev/null +++ b/Comp2D.R @@ -0,0 +1,143 @@ +## Time-stamp: "Last modified 2022-03-01 15:41:58 delucia" +library(ReacTran) +library(deSolve) +options(width=114) + +N <- 51 # number of grid cells +XX <- 1 # total size +dy <- dx <- XX/N # grid size +Dy <- Dx <- 0.1 # diffusion coeff, X- and Y-direction + +N2 <- ceiling(N/2) + +# The model equations + +Diff2D <- function (t, y, parms) { + CONC <- matrix(nrow = N, ncol = N, y) + # CONC[25, 25] <- 1 + dCONC <- tran.2D(CONC, D.x = Dx, D.y = Dy, dx = dx, dy = dy)$dC + + return (list(dCONC)) + +} + +# initial condition: 0 everywhere, except in central point +y <- matrix(nrow = N, ncol = N, data = 0) +y[N2, N2] <- 1 # initial concentration in the central point... + +as.numeric(y)[1301] + +# solve for 8 time units +times <- 0:5 +out <- ode.2D (y = y, func = Diff2D, t = times, parms = NULL, + dim = c(N,N), lrw=155412) + + +image(out, ask = FALSE, mfrow = c(2, 3), main = paste("time", times), asp=1) + +str(out) + +## adi <- data.matrix(data.table::fread("./bu2d/src/test2D", header=FALSE)) +## str(adi) + +adi <- data.table::fread("./bu2d/src/tt", header=FALSE) + +attributes(adi) <- attributes(out) + +madi1 <- matrix(adi[1,-1], 51,51,byrow = TRUE) +madi2 <- matrix(adi[2,-1], 51,51,byrow = TRUE) +madi6 <- matrix(adi[6,-1], 51,51,byrow = TRUE) +mout6 <- matrix(out[6,-1], 51,51,byrow = TRUE) + +madi6[1300] + + +madi6 <- matrix(unlist(adi[6,-1]), 501,501,byrow = TRUE) +class(madi6) +image(madi6, asp=1) +contour(madi6, asp=1) + +par(mfrow = c(1,3)) +image(madi1, asp=1) +image(madi2, asp=1) +image(madi5, asp=1) + +range(out[6, -1]) + +x11() + +plot(adi[6, -1], out[6, -1], pch=4) +abline(0,1) + +range(o6 <- out[6, -1]) +range(o2 <- out[2, -1]) + +image(madi6) + +contour(madi6) + +contour(mout6) + + +library(ReacTran) + +##################################################### +## FIRST SIMPLEST MODEL: Fick second law for diffusion +# d^2(C(x,t))/dx^2 = D. d^2(C(x,t))/dx^2 +# x: position +# t: time +# D: Diffusion coefficient +# C(x,t): Concentration at position x and time t + +# Model definition + +Fick_model=function(t=0,y,parms=NULL) +{ + + tran= tran.1D(C=y,D=1000,dx=grid)$dC + return(list(dC=tran)) +} + +# Grid definition + +C <- dnorm(seq(-10,10,.1)) +L <- 200 # length of the diffusion space +N <- length(C) # number of grid layers +grid <- setup.grid.1D(x.up = 0,L = L, N = N) + + +# Initial conditions + simulation + +Fick_solved_20 <- ode.1D(y = C, times=seq(0,20),func=Fick_model,nspec=1,method="lsoda") +Fick_solved_20000 <- ode.1D(y = C, times=seq(0,20000),func=Fick_model,nspec=1,method="lsoda") + + +# Plot Results + +plot(Fick_solved_20[1,2:dim(Fick_solved_20)[2]],type="l",xlab="x",ylab="Concentration",lwd=2,ylim=c(0,.5)) +par(new=T) +plot(Fick_solved_20[2,2:dim(Fick_solved_20)[2]],type="l",xlab="x",ylab="Concentration",lwd=2,ylim=c(0,.5),col="blue") +par(new=T) +plot(Fick_solved_20000[2,2:dim(Fick_solved_20000)[2]],type="l",xlab="x",ylab="Concentration",lwd=2,ylim=c(0,.5),col="red") + +Fick_solved_20[2,2:dim(Fick_solved_20)[2]]==Fick_solved_20000[2,2:dim(Fick_solved_20000)[2]] + +mass <- function(time, state, params) { + with(as.list(c(state, params)), { + H2O <- K3 * H * O + dH <- -K1 * H2O + dO <- -K2 * H2O + dH2O <- H2O + list(c(dH, dO, dH2O)) + }) +} + +params <- c(K1 = 2, + K2 = 1, + K3 = 0.005) +state <- c(H = 200, + O = 100, + H2O = 0) +time <- seq(0,5) +out <- ode(state, time, mass, params) +plot(out) diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt new file mode 100644 index 0000000..c92d35b --- /dev/null +++ b/app/CMakeLists.txt @@ -0,0 +1,8 @@ +add_executable(1D main_1D.cpp) +target_link_libraries(1D PUBLIC diffusion) + +add_executable(2D main_2D.cpp) +target_link_libraries(2D PUBLIC diffusion) + +add_executable(Comp2D main_2D_mdl.cpp) +target_link_libraries(Comp2D PUBLIC diffusion) diff --git a/src/main.cpp b/app/main_1D.cpp similarity index 69% rename from src/main.cpp rename to app/main_1D.cpp index 522c5c6..f7ec59b 100644 --- a/src/main.cpp +++ b/app/main_1D.cpp @@ -1,10 +1,14 @@ #include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET +#include "BoundaryCondition.hpp" #include // for copy, max +#include #include #include // for std #include // for vector using namespace std; +using namespace Diffusion; + int main(int argc, char *argv[]) { // dimension of grid @@ -15,6 +19,7 @@ int main(int argc, char *argv[]) { // create input + diffusion coefficients for each grid cell std::vector alpha(n, 1 * pow(10, -1)); std::vector field(n, 1 * std::pow(10, -6)); + std::vector bc(n, {0,0}); // create instance of diffusion module BTCSDiffusion diffu(dim); @@ -22,8 +27,9 @@ int main(int argc, char *argv[]) { diffu.setXDimensions(1, n); // set the boundary condition for the left ghost cell to dirichlet - diffu.setBoundaryCondition(0, BTCSDiffusion::BC_CONSTANT, - 5. * std::pow(10, -6)); + bc[0] = {Diffusion::BC_CONSTANT, 5*std::pow(10,-6)}; + // diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT, + // 5. * std::pow(10, -6)); // set timestep for simulation to 1 second diffu.setTimestep(1.); @@ -33,12 +39,12 @@ int main(int argc, char *argv[]) { // loop 100 times // output is currently generated by the method itself for (int i = 0; i < 100; i++) { - diffu.simulate(field, alpha); + diffu.simulate(field.data(), alpha.data(), bc.data()); cout << "Iteration: " << i << "\n\n"; - for (int j = 0; j < field.size(); j++) { - cout << field[j] << "\n"; + for (auto & cell : field) { + cout << cell << "\n"; } cout << "\n" << endl; diff --git a/app/main_2D.cpp b/app/main_2D.cpp new file mode 100644 index 0000000..e67004c --- /dev/null +++ b/app/main_2D.cpp @@ -0,0 +1,57 @@ +#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET +#include "BoundaryCondition.hpp" +#include // for copy, max +#include +#include +#include // for std +#include // for vector +using namespace std; + +using namespace Diffusion; + +int main(int argc, char *argv[]) { + + // dimension of grid + int dim = 2; + + int n = 5; + int m = 5; + + // create input + diffusion coefficients for each grid cell + std::vector alpha(n * m, 1 * pow(10, -1)); + std::vector field(n * m, 1 * std::pow(10, -6)); + std::vector bc(n*m, {0,0}); + + // create instance of diffusion module + BTCSDiffusion diffu(dim); + + diffu.setXDimensions(1, n); + diffu.setYDimensions(1, m); + + //set inlet to higher concentration without setting bc + field[12] = 5*std::pow(10,-3); + + // set timestep for simulation to 1 second + diffu.setTimestep(1.); + + cout << setprecision(12); + + for (int t = 0; t < 10; t++) { + diffu.simulate(field.data(), alpha.data(), bc.data()); + + cout << "Iteration: " << t << "\n\n"; + + // iterate through rows + for (int i = 0; i < m; i++) { + // iterate through columns + for (int j = 0; j < n; j++) { + cout << left << std::setw(20) << field[i * n + j]; + } + cout << "\n"; + } + + cout << "\n" << endl; + } + + return 0; +} diff --git a/app/main_2D_mdl.cpp b/app/main_2D_mdl.cpp new file mode 100644 index 0000000..d9b73ff --- /dev/null +++ b/app/main_2D_mdl.cpp @@ -0,0 +1,66 @@ +#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET +#include "BoundaryCondition.hpp" +#include // for copy, max +#include +#include +#include // for std +#include // for vector +using namespace std; + +using namespace Diffusion; + +int main(int argc, char *argv[]) { + + // dimension of grid + int dim = 2; + + int n = 501; + int m = 501; + + // create input + diffusion coefficients for each grid cell + std::vector alpha(n * m, 1 * pow(10, -1)); + std::vector field(n * m, 0.); + std::vector bc(n*m, {0,0}); + + field[125500] = 1; + + // create instance of diffusion module + BTCSDiffusion diffu(dim); + + diffu.setXDimensions(1., n); + diffu.setYDimensions(1., m); + + // set the boundary condition for the left ghost cell to dirichlet + //diffu.setBoundaryCondition(250, 250, BTCSDiffusion::BC_CONSTANT, 1); + // for (int d=0; d<5;d++){ + // diffu.setBoundaryCondition(d, 0, BC_CONSTANT, .1); + // } + // diffu.setBoundaryCondition(1, 1, BTCSDiffusion::BC_CONSTANT, .1); + // diffu.setBoundaryCondition(1, 1, BTCSDiffusion::BC_CONSTANT, .1); + + // set timestep for simulation to 1 second + diffu.setTimestep(1.); + + cout << setprecision(7); + + // First we output the initial state + cout << 0; + + for (int i=0; i < m*n; i++) { + cout << "," << field[i]; + } + cout << endl; + + // Now we simulate and output 8 steps à 1 sec + for (int t = 1; t < 6; t++) { + diffu.simulate(field.data(), alpha.data(), bc.data()); + cout << t; + + for (int i=0; i < m*n; i++) { + cout << "," << field[i]; + } + cout << endl; + } + + return 0; +} diff --git a/doc/.gitignore b/doc/.gitignore new file mode 100644 index 0000000..b51a70d --- /dev/null +++ b/doc/.gitignore @@ -0,0 +1,2 @@ +*.pdf +*.tex diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org new file mode 100644 index 0000000..6a86caf --- /dev/null +++ b/doc/ADI_scheme.org @@ -0,0 +1,106 @@ +#+TITLE: Adi 2D Scheme +#+LaTeX_CLASS_OPTIONS: [a4paper,10pt] +#+LATEX_HEADER: \usepackage{fullpage} +#+OPTIONS: toc:nil + +* Input + +- =c= $\rightarrow c$ + - containing current concentrations at each grid cell for species + - size: $N \times M$ + - row-major +- =alpha= $\rightarrow \alpha$ + - diffusion coefficient for both directions (x and y) + - size: $N \times M$ + - row-major +- =boundary_condition= $\rightarrow bc$ + - Defines closed or constant boundary condition for each grid cell + - size: $N \times M$ + - row-major + +* Internals + +- =A_matrix= $\rightarrow A$ + - coefficient matrix for linear equation system implemented as sparse matrix + - size: $((N+2)\cdot M) \times ((N+2)\cdot M)$ (including ghost zones in x direction) + - column-major (not relevant) + +- =b_vector= $\rightarrow b$ + - right hand side of the linear equation system + - size: $(N+2) \cdot M$ + - column-major (not relevant) +- =x_vector= $\rightarrow x$ + - solutions of the linear equation system + - size: $(N+2) \cdot M$ + - column-major (not relevant) + +* Calculation for $\frac{1}{2}$ timestep + +** Symbolic addressing of grid cells +[[./grid.png]] + +** Filling of matrix $A$ + +- row-wise iterating with $i$ over =c= and =\alpha= matrix respectively +- addressing each element of a row with $j$ +- matrix $A$ also containing $+2$ ghost nodes for each row of input matrix $\alpha$ + - $\rightarrow offset = N+2$ + - addressing each object $(i,j)$ in matrix $A$ with $(offset \cdot i + j, offset \cdot i + j)$ + +*** Rules + +$s_x(i,j) = \frac{\alpha(i,j)*\frac{t}{2}}{\Delta x^2}$ where $x$ defining the domain size in x direction. + +For the sake of simplicity we assume that each row of the $A$ matrix is addressed correctly with the given offset. + +**** Ghost nodes + +$A(i,-1) = 1$ + +$A(i,N) = 1$ + +**** Inlet + +$A(i,j) = \begin{cases} +1 & \text{if } bc(i,j) = \text{constant} \\ +-1-2*s_x(i,j) & \text{else} +\end{cases}$ + +$A(i,j\pm 1) = \begin{cases} +0 & \text{if } bc(i,j) = \text{constant} \\ +s_x(i,j) & \text{else} +\end{cases}$ + +** Filling of vector $b$ + +- each elements assign a concrete value to the according value of the row of matrix $A$ +- Adressing would look like this: $(i,j) = b(i \cdot (N+2) + j)$ + - $\rightarrow$ for simplicity we will write $b(i,j)$ + + + + +*** Rules + +**** Ghost nodes + +$b(i,-1) = \begin{cases} +0 & \text{if } bc(i,0) = \text{constant} \\ +c(i,0) & \text{else} +\end{cases}$ + +$b(i,N) = \begin{cases} +0 & \text{if } bc(i,N-1) = \text{constant} \\ +c(i,N-1) & \text{else} +\end{cases}$ + +*** Inlet + +$p(i,j) = \frac{\Delta t}{2}\alpha(i,j)\frac{c(i-1,j) - 2\cdot c(i,j) + c(i+1,j)}{\Delta x^2}$[fn:1] + +$b(i,j) = \begin{cases} +bc(i,j).\text{value} & \text{if } bc(i,N-1) = \text{constant} \\ +-c(i,j)-p(i,j) & \text{else} +\end{cases}$ + +[fn:1] $p$ is called =t0_c= inside code diff --git a/doc/grid.png b/doc/grid.png new file mode 100644 index 0000000..86940a4 Binary files /dev/null and b/doc/grid.png differ diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index da10b26..a9cd375 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -1,170 +1,324 @@ #include "BTCSDiffusion.hpp" +#include "BoundaryCondition.hpp" #include +#include +#include #include +#include #include +#include #include +#include +#include #include #include -const int BTCSDiffusion::BC_CONSTANT = 0; -const int BTCSDiffusion::BC_CLOSED = 1; -const int BTCSDiffusion::BC_FLUX = 2; +#include -BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) { - assert(dim <= 3); +constexpr int BTCS_MAX_DEP_PER_CELL = 3; +constexpr int BTCS_2D_DT_SIZE = 2; +constexpr double center_eq(double sx) { return -1. - 2. * sx; } + +Diffusion::BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) { grid_cells.resize(dim, 1); domain_size.resize(dim, 1); deltas.resize(dim, 1); + + this->time_step = 0; } -void BTCSDiffusion::setXDimensions(double domain_size, - unsigned int n_grid_cells) { - assert(this->grid_dim > 0); +void Diffusion::BTCSDiffusion::setXDimensions(double domain_size, + unsigned int n_grid_cells) { this->domain_size[0] = domain_size; this->grid_cells[0] = n_grid_cells; updateInternals(); } -void BTCSDiffusion::setYDimensions(double domain_size, - unsigned int n_grid_cells) { - assert(this->grid_dim > 1); +void Diffusion::BTCSDiffusion::setYDimensions(double domain_size, + unsigned int n_grid_cells) { this->domain_size[1] = domain_size; this->grid_cells[1] = n_grid_cells; updateInternals(); } -void BTCSDiffusion::setZDimensions(double domain_size, - unsigned int n_grid_cells) { - assert(this->grid_dim > 2); +void Diffusion::BTCSDiffusion::setZDimensions(double domain_size, + unsigned int n_grid_cells) { this->domain_size[2] = domain_size; this->grid_cells[2] = n_grid_cells; updateInternals(); } -unsigned int BTCSDiffusion::getXGridCellsN() { return this->grid_cells[0]; } -unsigned int BTCSDiffusion::getYGridCellsN() { return this->grid_cells[1]; } -unsigned int BTCSDiffusion::getZGridCellsN() { return this->grid_cells[2]; } -unsigned int BTCSDiffusion::getXDomainSize() { return this->domain_size[0]; } -unsigned int BTCSDiffusion::getYDomainSize() { return this->domain_size[1]; } -unsigned int BTCSDiffusion::getZDomainSize() { return this->domain_size[2]; } +auto Diffusion::BTCSDiffusion::getXGridCellsN() -> unsigned int { + return this->grid_cells[0]; +} +auto Diffusion::BTCSDiffusion::getYGridCellsN() -> unsigned int { + return this->grid_cells[1]; +} +auto Diffusion::BTCSDiffusion::getZGridCellsN() -> unsigned int { + return this->grid_cells[2]; +} +auto Diffusion::BTCSDiffusion::getXDomainSize() -> double { + return this->domain_size[0]; +} +auto Diffusion::BTCSDiffusion::getYDomainSize() -> double { + return this->domain_size[1]; +} +auto Diffusion::BTCSDiffusion::getZDomainSize() -> double { + return this->domain_size[2]; +} -void BTCSDiffusion::updateInternals() { +void Diffusion::BTCSDiffusion::updateInternals() { for (int i = 0; i < grid_dim; i++) { deltas[i] = (double)domain_size[i] / grid_cells[i]; } +} +void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c, + const BCVectorRowMajor &bc, + const DVectorRowMajor &alpha, + double dx, double time_step, + int size, + const DVectorRowMajor &t0_c) { - int cells = 1; + reserveMemory(size, BTCS_MAX_DEP_PER_CELL); - for (int i = 0; i < grid_dim; i++) { - cells *= (grid_cells[i] + 2); - } + fillMatrixFromRow(alpha.row(0), bc.row(0), size, dx, time_step); + fillVectorFromRow(c, alpha, bc, Eigen::VectorXd::Constant(size, 0), size, dx, + time_step); - bc.resize(cells, {BTCSDiffusion::BC_CLOSED, 0}); + solveLES(); + + c = x_vector.segment(1, size); } -void BTCSDiffusion::simulate1D(std::vector &c, boundary_condition left, - boundary_condition right, - const std::vector &alpha, double dx, - int size) { +inline void Diffusion::BTCSDiffusion::reserveMemory(int size, + int max_count_per_line) { + size += 2; - bool left_is_constant = (left.type == BTCSDiffusion::BC_CONSTANT); - bool right_is_constant = (right.type == BTCSDiffusion::BC_CONSTANT); + A_matrix.resize(size, size); + A_matrix.reserve(Eigen::VectorXi::Constant(size, max_count_per_line)); - //The sizes for matrix and vectors of the equation system is defined by the - //actual size of the input vector and if the system is (partially) closed. - //Then we will need ghost nodes. So this variable will give the count of ghost - //nodes. - int bc_offset = !left_is_constant + !right_is_constant; - ; + b_vector.resize(size); + x_vector.resize(size); +} - // set sizes of private and yet allocated vectors - b_vector.resize(size + bc_offset); - x_vector.resize(size + bc_offset); +void Diffusion::BTCSDiffusion::simulate1D( + Eigen::Map &c, Eigen::Map &alpha, + Eigen::Map &bc) { - /* - * Begin to solve the equation system using LU solver of Eigen. - * - * But first fill the A matrix and b vector. - */ + int size = this->grid_cells[0]; + double dx = this->deltas[0]; + double time_step = this->time_step; - // Set boundary condition for ghost nodes (for closed or flux system) or outer - // inlet nodes (constant boundary condition) - A_matrix.resize(size + bc_offset, size + bc_offset); - A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3)); + DVectorRowMajor input_field = c.row(0); + + simulate_base(input_field, bc, alpha, dx, time_step, size, + Eigen::VectorXd::Constant(size, 0)); + + c.row(0) << input_field; +} + +void Diffusion::BTCSDiffusion::simulate2D( + Eigen::Map &c, Eigen::Map &alpha, + Eigen::Map &bc) { + + int n_rows = this->grid_cells[1]; + int n_cols = this->grid_cells[0]; + double dx = this->deltas[0]; + DMatrixRowMajor t0_c; + + double local_dt = this->time_step / BTCS_2D_DT_SIZE; + + t0_c = calc_t0_c(c, alpha, bc, local_dt, dx); + + for (int i = 0; i < n_rows; i++) { + DVectorRowMajor input_field = c.row(i); + simulate_base(input_field, bc.row(i), alpha.row(i), dx, local_dt, n_cols, + t0_c.row(i)); + c.row(i) << input_field; + } + + dx = this->deltas[1]; + + t0_c = + calc_t0_c(c.transpose(), alpha.transpose(), bc.transpose(), local_dt, dx); + + for (int i = 0; i < n_cols; i++) { + DVectorRowMajor input_field = c.col(i); + simulate_base(input_field, bc.col(i), alpha.col(i), dx, local_dt, n_rows, + t0_c.row(i)); + c.col(i) << input_field.transpose(); + } +} + +auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c, + const DMatrixRowMajor &alpha, + const BCMatrixRowMajor &bc, + double time_step, double dx) + -> DMatrixRowMajor { + + int n_rows = this->grid_cells[1]; + int n_cols = this->grid_cells[0]; + + DMatrixRowMajor t0_c(n_rows, n_cols); + + std::array y_values; + + // first, iterate over first row + for (int j = 0; j < n_cols; j++) { + y_values[0] = getBCFromFlux(bc(0, j), c(0, j), alpha(0, j)); + y_values[1] = c(0, j); + y_values[2] = c(1, j); + + t0_c(0, j) = time_step * alpha(0, j) * + (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); + } + + // then iterate over inlet + for (int i = 1; i < n_rows - 1; i++) { + for (int j = 0; j < n_cols; j++) { + + y_values[0] = c(i - 1, j); + y_values[1] = c(i, j); + y_values[2] = c(i + 1, j); + + t0_c(i, j) = time_step * alpha(i, j) * + (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); + } + } + + int end = n_rows - 1; + + // and finally over last row + for (int j = 0; j < n_cols; j++) { + y_values[0] = c(end - 1, j); + y_values[1] = c(end, j); + y_values[2] = getBCFromFlux(bc(end, j), c(end, j), alpha(end, j)); + + t0_c(end, j) = time_step * alpha(end, j) * + (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); + } + + return t0_c; +} + +inline void Diffusion::BTCSDiffusion::fillMatrixFromRow( + const DVectorRowMajor &alpha, const BCVectorRowMajor &bc, int size, + double dx, double time_step) { + + Diffusion::boundary_condition left = bc[0]; + Diffusion::boundary_condition right = bc[size - 1]; + + bool left_constant = (left.type == Diffusion::BC_CONSTANT); + bool right_constant = (right.type == Diffusion::BC_CONSTANT); + + int A_size = A_matrix.cols(); A_matrix.insert(0, 0) = 1; - b_vector[0] = - (left_is_constant ? left.value : getBCFromFlux(left, c[0], alpha[0])); - A_matrix.insert((size + bc_offset) - 1, (size + bc_offset) - 1) = 1; - b_vector[size + bc_offset - 1] = - (right_is_constant ? right.value - : getBCFromFlux(right, c[size - 1], alpha[size - 1])); + if (left_constant) { + A_matrix.insert(1, 1) = 1; + } - // Start filling the A matrix - // =i= is used for equation system matrix and vector indexing - // and =j= for indexing of c,alpha and bc - for (int i = 1, j = i + !(left_is_constant); i < size - right_is_constant; - i++, j++) { + A_matrix.insert(A_size - 1, A_size - 1) = 1; - // if current grid cell is considered as constant boundary conditon - if (bc[j].type == BTCSDiffusion::BC_CONSTANT) { - A_matrix.insert(i, i) = 1; - b_vector[i] = bc[j].value; + if (right_constant) { + A_matrix.insert(A_size - 2, A_size - 2) = 1; + } + + for (int j = 1 + (int)left_constant, k = j - 1; + k < size - (int)right_constant; j++, k++) { + double sx = (alpha[k] * time_step) / (dx * dx); + + if (bc[k].type == Diffusion::BC_CONSTANT) { + A_matrix.insert(j, j) = 1; continue; } - double sx = (alpha[j] * time_step) / (dx * dx); - - A_matrix.insert(i, i) = -1. - 2. * sx; - A_matrix.insert(i, i - 1) = sx; - A_matrix.insert(i, i + 1) = sx; - - b_vector[i] = -c[j]; - } - - // start to solve - Eigen::SparseLU, Eigen::COLAMDOrdering> - solver; - solver.analyzePattern(A_matrix); - - solver.factorize(A_matrix); - - x_vector = solver.solve(b_vector); - - //fill solution back in place into =c= vector - for (int i = 0, j = i + !left_is_constant; i < c.size(); i++, j++) { - c[i] = x_vector[i + !left_is_constant]; + A_matrix.insert(j, j) = center_eq(sx); + A_matrix.insert(j, (j - 1)) = sx; + A_matrix.insert(j, (j + 1)) = sx; } } -void BTCSDiffusion::setTimestep(double time_step) { +inline void Diffusion::BTCSDiffusion::fillVectorFromRow( + const DVectorRowMajor &c, const DVectorRowMajor &alpha, + const BCVectorRowMajor &bc, const DVectorRowMajor &t0_c, int size, + double dx, double time_step) { + + Diffusion::boundary_condition left = bc[0]; + Diffusion::boundary_condition right = bc[size - 1]; + + bool left_constant = (left.type == Diffusion::BC_CONSTANT); + bool right_constant = (right.type == Diffusion::BC_CONSTANT); + + int b_size = b_vector.size(); + + for (int j = 0; j < size; j++) { + boundary_condition tmp_bc = bc[j]; + + if (tmp_bc.type == Diffusion::BC_CONSTANT) { + b_vector[j + 1] = tmp_bc.value; + continue; + } + + double t0_c_j = time_step * alpha[j] * (t0_c[j] / (dx * dx)); + b_vector[j + 1] = -c[j] - t0_c_j; + } + + if (!left_constant) { + // this is not correct currently.We will fix this when we are able to define + // FLUX boundary conditions + b_vector[0] = getBCFromFlux(left, c[0], alpha[0]); + } + + if (!right_constant) { + b_vector[b_size - 1] = getBCFromFlux(right, c[size - 1], alpha[size - 1]); + } +} + +void Diffusion::BTCSDiffusion::setTimestep(double time_step) { this->time_step = time_step; } -void BTCSDiffusion::simulate(std::vector &c, - const std::vector &alpha) { +void Diffusion::BTCSDiffusion::simulate(double *c, double *alpha, + Diffusion::boundary_condition *bc) { if (this->grid_dim == 1) { - simulate1D(c, bc[0], bc[grid_cells[0] + 1], alpha, this->deltas[0], - this->grid_cells[0]); + Eigen::Map c_in(c, this->grid_cells[0]); + Eigen::Map alpha_in(alpha, this->grid_cells[0]); + Eigen::Map bc_in(bc, this->grid_cells[0]); + + simulate1D(c_in, alpha_in, bc_in); + } + if (this->grid_dim == 2) { + Eigen::Map c_in(c, this->grid_cells[1], + this->grid_cells[0]); + + Eigen::Map alpha_in(alpha, this->grid_cells[1], + this->grid_cells[0]); + + Eigen::Map bc_in(bc, this->grid_cells[1], + this->grid_cells[0]); + + simulate2D(c_in, alpha_in, bc_in); } } -inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc, - double neighbor_c, - double neighbor_alpha) { +inline auto Diffusion::BTCSDiffusion::getBCFromFlux(boundary_condition bc, + double neighbor_c, + double neighbor_alpha) + -> double { - double val; + double val = 0; - if (bc.type == BTCSDiffusion::BC_CLOSED) { + if (bc.type == Diffusion::BC_CLOSED) { val = neighbor_c; - } else if (bc.type == BTCSDiffusion::BC_FLUX) { + } else if (bc.type == Diffusion::BC_FLUX) { // TODO // val = bc[index].value; } else { @@ -174,8 +328,13 @@ inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc, return val; } -void BTCSDiffusion::setBoundaryCondition(int index, bctype type, double value) { +inline void Diffusion::BTCSDiffusion::solveLES() { + // start to solve + Eigen::SparseLU, Eigen::COLAMDOrdering> + solver; + solver.analyzePattern(A_matrix); - bc[index].type = type; - bc[index].value = value; + solver.factorize(A_matrix); + + x_vector = solver.solve(b_vector); } diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 7614ce5..b994b80 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -1,16 +1,18 @@ #ifndef BTCSDIFFUSION_H_ #define BTCSDIFFUSION_H_ +#include "BoundaryCondition.hpp" + #include +#include +#include +#include +#include #include #include #include -/*! - * Defines both types of boundary condition as a datatype. - */ -typedef int bctype; - +namespace Diffusion { /*! * Class implementing a solution for a 1/2/3D diffusion equation using backward * euler. @@ -18,21 +20,6 @@ typedef int bctype; class BTCSDiffusion { public: - /*! - * Defines a constant/Dirichlet boundary condition. - */ - static const int BC_CONSTANT; - - /*! - * Defines a closed/Neumann boundary condition. - */ - static const int BC_CLOSED; - - /*! - * Defines a flux/Cauchy boundary condition. - */ - static const int BC_FLUX; - /*! * Creates a diffusion module. * @@ -75,37 +62,37 @@ public: /*! * Returns the number of grid cells in x direction. */ - unsigned int getXGridCellsN(); + auto getXGridCellsN() -> unsigned int; /*! * Returns the number of grid cells in y direction. */ - unsigned int getYGridCellsN(); + auto getYGridCellsN() -> unsigned int; /*! * Returns the number of grid cells in z direction. */ - unsigned int getZGridCellsN(); + auto getZGridCellsN() -> unsigned int; /*! * Returns the domain size in x direction. */ - unsigned int getXDomainSize(); + auto getXDomainSize() -> double; /*! * Returns the domain size in y direction. */ - unsigned int getYDomainSize(); + auto getYDomainSize() -> double; /*! * Returns the domain size in z direction. */ - unsigned int getZDomainSize(); + auto getZDomainSize() -> double; /*! * With given ghost zones simulate diffusion. Only 1D allowed at this moment. * - * @param c Vector describing the concentration of one solution of the grid as - * continious memory (row major). - * @param alpha Vector of diffusion coefficients for each grid element. + * @param c Pointer to continious memory describing the current concentration state of each grid cell. + * @param alpha Pointer to memory area of diffusion coefficients for each grid element. + * @param bc Pointer to memory setting boundary conidition of each grid cell. */ - void simulate(std::vector &c, const std::vector &alpha); + void simulate(double *c, double *alpha, Diffusion::boundary_condition *bc); /*! * Set the timestep of the simulation @@ -114,48 +101,62 @@ public: */ void setTimestep(double time_step); - /*! - * Set the boundary condition of the given grid. This is done by defining an - * index (exact order still to be determined), the type of the boundary - * condition and the according value. - * - * @param index Index of the grid cell the boundary condition is applied to. - * @param type Type of the boundary condition. Must be constant, closed or - * flux. - * @param value For constant boundary conditions this value is set - * during solving. For flux value refers to a gradient of change for this grid - * cell. For closed this value has no effect since a gradient of 0 is used. - */ - void setBoundaryCondition(int index, bctype type, double value); - private: - typedef struct boundary_condition { - bctype type; - double value; - } boundary_condition; - typedef Eigen::Triplet T; + typedef Eigen::Matrix + DMatrixRowMajor; + typedef Eigen::Matrix + DVectorRowMajor; + typedef Eigen::Matrix + BCMatrixRowMajor; + typedef Eigen::Matrix + BCVectorRowMajor; - void simulate1D(std::vector &c, boundary_condition left, - boundary_condition right, const std::vector &alpha, - double dx, int size); - void simulate2D(std::vector &c); + void simulate_base(DVectorRowMajor &c, const BCVectorRowMajor &bc, + const DVectorRowMajor &alpha, double dx, double time_step, + int size, const DVectorRowMajor &t0_c); + + void simulate1D(Eigen::Map &c, + Eigen::Map &alpha, + Eigen::Map &bc); + + void simulate2D(Eigen::Map &c, + Eigen::Map &alpha, + Eigen::Map &bc); + + auto calc_t0_c(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha, + const BCMatrixRowMajor &bc, double time_step, double dx) + -> DMatrixRowMajor; + + inline void fillMatrixFromRow(const DVectorRowMajor &alpha, + const BCVectorRowMajor &bc, int size, double dx, + double time_step); + inline void fillVectorFromRow(const DVectorRowMajor &c, + const DVectorRowMajor &alpha, + const BCVectorRowMajor &bc, + const DVectorRowMajor &t0_c, int size, + double dx, double time_step); void simulate3D(std::vector &c); - inline double getBCFromFlux(boundary_condition bc, double nearest_value, - double neighbor_alpha); - void updateInternals(); - std::vector bc; + inline void reserveMemory(int size, int max_count_per_line); + inline static auto getBCFromFlux(Diffusion::boundary_condition bc, + double neighbor_c, double neighbor_alpha) + -> double; + + void solveLES(); + void updateInternals(); Eigen::SparseMatrix A_matrix; Eigen::VectorXd b_vector; Eigen::VectorXd x_vector; double time_step; - int grid_dim; + unsigned int grid_dim; std::vector grid_cells; - std::vector domain_size; + std::vector domain_size; std::vector deltas; }; - +} // namespace Diffusion #endif // BTCSDIFFUSION_H_ diff --git a/src/BoundaryCondition.hpp b/src/BoundaryCondition.hpp new file mode 100644 index 0000000..e47482d --- /dev/null +++ b/src/BoundaryCondition.hpp @@ -0,0 +1,33 @@ +#ifndef BOUNDARYCONDITION_H_ +#define BOUNDARYCONDITION_H_ + +namespace Diffusion { + +/*! + * Defines both types of boundary condition as a datatype. + */ +typedef int bctype; + +/*! + * Defines a closed/Neumann boundary condition. + */ +static const bctype BC_CLOSED = 0; + +/*! + * Defines a flux/Cauchy boundary condition. + */ +static const bctype BC_FLUX = 1; + +/*! + * Defines a constant/Dirichlet boundary condition. + */ +static const bctype BC_CONSTANT = 2; + +typedef struct boundary_condition { + bctype type; + double value; +} boundary_condition; + +} // namespace Diffusion + +#endif // BOUNDARYCONDITION_H_ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fab0508..ffc3958 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,5 +1,3 @@ add_library(diffusion OBJECT BTCSDiffusion.cpp BTCSDiffusion.hpp) target_link_libraries(diffusion Eigen3::Eigen) - -add_executable(test main.cpp) -target_link_libraries(test PUBLIC diffusion) +target_include_directories(diffusion PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})