max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-13 17:38:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch '2D' into 'main'

Browse Source 
Implementing 2D ADI-BTCS

See merge request mluebke/diffusion!5
This commit is contained in:
Max Lübke 2022-03-08 15:05:22 +01:00
commit a24c26beb0
15 changed files with 1011 additions and 170 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
.gitlab-ci.yml
View File

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

1
CMakeLists.txt
View File

@ -7,4 +7,5 @@ set(CMAKE_CXX_STANDARD 14)
find_package(Eigen3 REQUIRED NO_MODULE)
add_subdirectory(app)
add_subdirectory(src)

227
Comp.R Normal file
View File

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

143
Comp2D.R Normal file
View File

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

8
app/CMakeLists.txt Normal file
View File

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

16
src/main.cpp → app/main_1D.cpp
View File

@ -1,10 +1,14 @@
#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
#include "BoundaryCondition.hpp"
#include <algorithm> // for copy, max
#include <cmath>
#include <iomanip>
#include <iostream> // for std
#include <vector> // 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<double> alpha(n, 1 * pow(10, -1));
std::vector<double> field(n, 1 * std::pow(10, -6));
std::vector<boundary_condition> 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;

57
app/main_2D.cpp Normal file
View File

@ -0,0 +1,57 @@
#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
#include "BoundaryCondition.hpp"
#include <algorithm> // for copy, max
#include <cmath>
#include <iomanip>
#include <iostream> // for std
#include <vector> // 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<double> alpha(n * m, 1 * pow(10, -1));
std::vector<double> field(n * m, 1 * std::pow(10, -6));
std::vector<boundary_condition> 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;
}

66
app/main_2D_mdl.cpp Normal file
View File

@ -0,0 +1,66 @@
#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
#include "BoundaryCondition.hpp"
#include <algorithm> // for copy, max
#include <cmath>
#include <iomanip>
#include <iostream> // for std
#include <vector> // 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<double> alpha(n * m, 1 * pow(10, -1));
std::vector<double> field(n * m, 0.);
std::vector<boundary_condition> 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;
}

2
doc/.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
*.pdf
*.tex

106
doc/ADI_scheme.org Normal file
View File

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

BIN
doc/grid.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 151 KiB

347
src/BTCSDiffusion.cpp
View File

@ -1,170 +1,324 @@
#include "BTCSDiffusion.hpp"
#include "BoundaryCondition.hpp"
#include <Eigen/SparseLU>
#include <Eigen/src/Core/Map.h>
#include <Eigen/src/Core/Matrix.h>
#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <iomanip>
#include <iterator>
#include <ostream>
#include <tuple>
#include <vector>
const int BTCSDiffusion::BC_CONSTANT = 0;
const int BTCSDiffusion::BC_CLOSED = 1;
const int BTCSDiffusion::BC_FLUX = 2;
#include <iostream>
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,
void Diffusion::BTCSDiffusion::setXDimensions(double domain_size,
unsigned int n_grid_cells) {
assert(this->grid_dim > 0);
this->domain_size[0] = domain_size;
this->grid_cells[0] = n_grid_cells;
updateInternals();
}
void BTCSDiffusion::setYDimensions(double domain_size,
void Diffusion::BTCSDiffusion::setYDimensions(double domain_size,
unsigned int n_grid_cells) {
assert(this->grid_dim > 1);
this->domain_size[1] = domain_size;
this->grid_cells[1] = n_grid_cells;
updateInternals();
}
void BTCSDiffusion::setZDimensions(double domain_size,
void Diffusion::BTCSDiffusion::setZDimensions(double domain_size,
unsigned int n_grid_cells) {
assert(this->grid_dim > 2);
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);
solveLES();
c = x_vector.segment(1, size);
}
bc.resize(cells, {BTCSDiffusion::BC_CLOSED, 0});
inline void Diffusion::BTCSDiffusion::reserveMemory(int size,
int max_count_per_line) {
size += 2;
A_matrix.resize(size, size);
A_matrix.reserve(Eigen::VectorXi::Constant(size, max_count_per_line));
b_vector.resize(size);
x_vector.resize(size);
}
void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
boundary_condition right,
const std::vector<double> &alpha, double dx,
int size) {
void Diffusion::BTCSDiffusion::simulate1D(
Eigen::Map<DVectorRowMajor> &c, Eigen::Map<const DVectorRowMajor> &alpha,
Eigen::Map<const BCVectorRowMajor> &bc) {
bool left_is_constant = (left.type == BTCSDiffusion::BC_CONSTANT);
bool right_is_constant = (right.type == BTCSDiffusion::BC_CONSTANT);
int size = this->grid_cells[0];
double dx = this->deltas[0];
double time_step = this->time_step;
//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;
;
DVectorRowMajor input_field = c.row(0);
// set sizes of private and yet allocated vectors
b_vector.resize(size + bc_offset);
x_vector.resize(size + bc_offset);
simulate_base(input_field, bc, alpha, dx, time_step, size,
Eigen::VectorXd::Constant(size, 0));
/*
* Begin to solve the equation system using LU solver of Eigen.
*
* But first fill the A matrix and b vector.
*/
c.row(0) << input_field;
}
// 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));
void Diffusion::BTCSDiffusion::simulate2D(
Eigen::Map<DMatrixRowMajor> &c, Eigen::Map<const DMatrixRowMajor> &alpha,
Eigen::Map<const BCMatrixRowMajor> &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<double, 3> 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::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
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<double> &c,
const std::vector<double> &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],
Eigen::Map<DVectorRowMajor> c_in(c, this->grid_cells[0]);
Eigen::Map<const DVectorRowMajor> alpha_in(alpha, this->grid_cells[0]);
Eigen::Map<const BCVectorRowMajor> bc_in(bc, this->grid_cells[0]);
simulate1D(c_in, alpha_in, bc_in);
}
if (this->grid_dim == 2) {
Eigen::Map<DMatrixRowMajor> c_in(c, this->grid_cells[1],
this->grid_cells[0]);
Eigen::Map<const DMatrixRowMajor> alpha_in(alpha, this->grid_cells[1],
this->grid_cells[0]);
Eigen::Map<const BCMatrixRowMajor> 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,
inline auto Diffusion::BTCSDiffusion::getBCFromFlux(boundary_condition bc,
double neighbor_c,
double neighbor_alpha) {
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::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
solver;
solver.analyzePattern(A_matrix);
bc[index].type = type;
bc[index].value = value;
solver.factorize(A_matrix);
x_vector = solver.solve(b_vector);
}

121
src/BTCSDiffusion.hpp
View File

@ -1,16 +1,18 @@
#ifndef BTCSDIFFUSION_H_
#define BTCSDIFFUSION_H_
#include "BoundaryCondition.hpp"
#include <Eigen/Sparse>
#include <Eigen/src/Core/Map.h>
#include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/Core/util/Constants.h>
#include <cstddef>
#include <tuple>
#include <type_traits>
#include <vector>
/*!
* 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<double> &c, const std::vector<double> &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<double> T;
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
DMatrixRowMajor;
typedef Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>
DVectorRowMajor;
typedef Eigen::Matrix<Diffusion::boundary_condition, Eigen::Dynamic,
Eigen::Dynamic, Eigen::RowMajor>
BCMatrixRowMajor;
typedef Eigen::Matrix<Diffusion::boundary_condition, 1, Eigen::Dynamic,
Eigen::RowMajor>
BCVectorRowMajor;
void simulate1D(std::vector<double> &c, boundary_condition left,
boundary_condition right, const std::vector<double> &alpha,
double dx, int size);
void simulate2D(std::vector<double> &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<DVectorRowMajor> &c,
Eigen::Map<const DVectorRowMajor> &alpha,
Eigen::Map<const BCVectorRowMajor> &bc);
void simulate2D(Eigen::Map<DMatrixRowMajor> &c,
Eigen::Map<const DMatrixRowMajor> &alpha,
Eigen::Map<const BCMatrixRowMajor> &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<double> &c);
inline double getBCFromFlux(boundary_condition bc, double nearest_value,
double neighbor_alpha);
void updateInternals();
std::vector<boundary_condition> 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<double> A_matrix;
Eigen::VectorXd b_vector;
Eigen::VectorXd x_vector;
double time_step;
int grid_dim;
unsigned int grid_dim;
std::vector<unsigned int> grid_cells;
std::vector<unsigned int> domain_size;
std::vector<double> domain_size;
std::vector<double> deltas;
};
} // namespace Diffusion
#endif // BTCSDIFFUSION_H_

33
src/BoundaryCondition.hpp Normal file
View File

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

4
src/CMakeLists.txt
View File

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