diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 205a2e3..d9f4baf 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,4 +1,4 @@ -image: sobc/gitlab-ci +image: git.gfz-potsdam.de:5000/naaice/tug:ci stages: - build @@ -6,8 +6,6 @@ stages: - static_analyze build_release: - before_script: - - apt-get update && apt-get install -y libeigen3-dev git stage: build artifacts: paths: @@ -24,11 +22,11 @@ test: - ./build/test/testTug lint: - before_script: - - apt-get update && apt-get install -y libeigen3-dev stage: static_analyze + before_script: + - apk add clang-extra-tools script: - mkdir lint && cd lint - - cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*,readability-*, modernize-*" -DTUG_ENABLE_TESTING=OFF .. + - cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*, modernize-*" -DTUG_ENABLE_TESTING=OFF .. - make tug when: manual diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index e0d2abd..1124a92 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -1,7 +1,7 @@ -#+TITLE: Numerical solution of diffusion equation in 2D with ADI Scheme +#+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{amsmath, systeme} +#+LATEX_HEADER: \usepackage{amsmath, systeme, cancel, xcolor} #+OPTIONS: toc:nil @@ -308,8 +308,8 @@ form: * Heterogeneous diffusion -If the diffusion coefficient $\alpha$ is spatially variable, [[eqn:1]] can -be rewritten: +If the diffusion coefficient $\alpha$ is spatially variable, equation +[[eqn:1]] can be rewritten as: #+NAME: eqn:hetdiff \begin{align} @@ -391,8 +391,8 @@ concentrations. ** Direct discretization As noted in literature (LeVeque and Numerical Recipes) a better way is -to discretize directly the physical problem ([[eqn:hetdiff]]) at points -halfway between grid points: +to discretize directly the physical problem (eq. [[eqn:hetdiff]]) at +points halfway between grid points: \begin{align*} \begin{cases} @@ -401,15 +401,18 @@ halfway between grid points: \end{cases} \end{align*} -\noindent A further differentiation gives us the centered +\noindent A further differentiation gives us the spatially centered approximation of $\frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right)$: -\begin{align*} +#+NAME: eqn:CS_het +\begin{equation} +\begin{aligned} \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right)(x_i) & \simeq \frac{1}{\Delta x}\left[\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) - \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \right]\\ &\displaystyle =\frac{1}{\Delta x^2} \left[ \alpha_{i+1/2}C_{i+1} - (\alpha_{i+1/2}+\alpha_{i-1/2}) C_{i} + \alpha_{i-1/2}C_{i-1}\right] -\end{align*} +\end{aligned} +\end{equation} \noindent The ADI scheme with this approach becomes: @@ -439,3 +442,98 @@ explicit terms, the two sweeps become: \end{aligned} \right. \end{equation} + +\bigskip + +\noindent The "interblock" diffusion coefficients $\alpha_{i+1/2,j}$ +can be arithmetic mean: + +\[ +\displaystyle \alpha_{i+1/2, j} = \displaystyle \frac{\alpha_{i+1, j} + \alpha_{i, j}}{2} +\] + +\noindent or the harmonic mean: + +\[ +\displaystyle \alpha_{i+1/2, j} = \displaystyle \frac{2}{\frac{1}{\alpha_{i+1, j}} + \frac{1}{\alpha_{i, j}}} +\] + + + +\pagebreak + +* Explicit scheme for 2D heterogeneous diffusion + +A classical explicit FTCS scheme (forward in time, central in space) +for 2D heterogeneous diffusion can be expressed simply leveraging the +discretization of equation [[eqn:CS_het]]: + +#+NAME: eqn:2DHeterFTCS +\begin{equation} +\begin{aligned} +\frac{C_{i,j}^{t+1} - C_{i,j}^{t}}{\Delta t} = & \frac{1}{\Delta x^2} \left[ \alpha^x_{i+1/2, j}C^t_{i+1, j} - (\alpha^x_{i+1/2, j} + \alpha^x_{i-1/2, j}) C^t_{i,j} + \alpha^x_{i-1/2,j}C^t_{i-1,j}\right] + \\ + & \frac{1}{\Delta y^2} \left[ \alpha^y_{i, j+1/2}C^t_{i, j+1} - (\alpha^y_{i, j+1/2} + \alpha^y_{i, j-1/2}) C^t_{i,j} + \alpha^y_{i,j-1/2}C^t_{i,j-1}\right] +\end{aligned} +\end{equation} +\noindent where in the RHS only the known concentrations at time $t$ +appear. Rearranging the terms, we get: + +#+NAME: eqn:2DHeterFTCS_final +\begin{equation} +\begin{aligned} +C_{i,j}^{t+1} = & C^t_{i,j} +\\ + & \frac{\Delta t}{\Delta x^2} \left[ \alpha^x_{i+1/2, j}C^t_{i+1, j} - (\alpha^x_{i+1/2, j} + \alpha^x_{i-1/2, j}) C^t_{i,j} + \alpha^x_{i-1/2,j}C^t_{i-1,j}\right] + \\ + & \frac{\Delta t}{\Delta y^2} \left[ \alpha^y_{i, j+1/2}C^t_{i, j+1} - (\alpha^y_{i, j+1/2} + \alpha^y_{i, j-1/2}) C^t_{i,j} + \alpha^y_{i,j-1/2}C^t_{i,j-1}\right] +\end{aligned} +\end{equation} + +The Courant-Friedrichs-Lewy stability criterion for this scheme reads: +#+NAME: eqn:CFL2DFTCS +\begin{equation} +\Delta t \leq \frac{(\Delta x^2, \Delta y^2)}{2 \max(\alpha_{i,j})} +\end{equation} + +** Boundary conditions + +In analogy to the treatment of the 1D homogeneous FTCS scheme (cfr +section 1), we need to differentiate the domain boundaries ($i=0$ and +$i=n_x$; the same applies to $j$ of course) accounting for the +discrepancy in the discretization. + +For the zero-th (left) cell, whose center is at $x=dx/2$, we can +evaluate the left gradient with the left boundary using such distance, +calling $l$ the numerical value of a constant boundary condition, +equation [[eqn:CS_het]] becomes: + +#+NAME: eqn:2D_FTCS_left +\begin{equation} +\begin{aligned} +\frac{\partial}{\partial x} \left(\alpha(x) +\frac{\partial C }{\partial x}\right)(x_0) & \simeq \frac{1}{\Delta x}\left[\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) - \alpha_{i} \left( \frac{C_{i} - l }{\frac{\Delta x}{2}} \right) \right]\\ +&\displaystyle =\frac{1}{\Delta x^2} \left[ \alpha_{i+1/2}C_{i+1} - (\alpha_{i+1/2}+ 2\alpha_i) C_{i} + 2 \alpha_{i}\cdot l\right] +\end{aligned} +\end{equation} + +\noindent Similarly, for $i=n_x$, + +#+NAME: eqn:2D_FTCS_right +\begin{equation} +\begin{aligned} +\frac{\partial}{\partial x} \left(\alpha(x) +\frac{\partial C }{\partial x}\right)(x_n) & \simeq \frac{1}{\Delta x}\left[\alpha_{i} \left( \frac{r -C_{i}}{\frac{\Delta x}{2}} \right) - \alpha_{i-1/2} \left( \frac{C_{i} - C_{i-1}} {\Delta x} \right) \right]\\ +&\displaystyle =\frac{1}{\Delta x^2} \left[ 2 \alpha_{i} r - (\alpha_{i+1/2}+ 2\alpha_i) C_{i} + \alpha_{i-1/2}\cdot C_{i-1}\right] +\end{aligned} +\end{equation} + +If on the right boundary we have *closed* or Neumann condition, the +left derivative becomes zero and we are left with: + +#+NAME: eqn:2D_FTCS_rightclosed +\begin{equation} +\begin{aligned} +\frac{\partial}{\partial x} \left(\alpha(x) +\frac{\partial C }{\partial x}\right)(x_n) & \simeq \frac{1}{\Delta x}\left[\cancel{\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right)} - \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \right]\\ +&\displaystyle =\frac{\alpha_{i-1/2}}{\Delta x^2} (C_{i-1} - C_i) +\end{aligned} +\end{equation} + diff --git a/doc/ADI_scheme.pdf b/doc/ADI_scheme.pdf index 41b0f4d..7a2fddb 100644 Binary files a/doc/ADI_scheme.pdf and b/doc/ADI_scheme.pdf differ diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index b806fac..47042c2 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -2,9 +2,11 @@ add_executable(first_example first_example.cpp) add_executable(second_example second_example.cpp) add_executable(boundary_example1D boundary_example1D.cpp) add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp) +add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp) target_link_libraries(first_example tug) target_link_libraries(second_example tug) target_link_libraries(boundary_example1D tug) target_link_libraries(FTCS_2D_proto_example tug) +target_link_libraries(FTCS_1D_proto_example tug) target_link_libraries(FTCS_2D_proto_example easy_profiler) \ No newline at end of file diff --git a/examples/FTCS_1D_proto_example.cpp b/examples/FTCS_1D_proto_example.cpp new file mode 100644 index 0000000..c0936c4 --- /dev/null +++ b/examples/FTCS_1D_proto_example.cpp @@ -0,0 +1,44 @@ +#include + +int main(int argc, char *argv[]) { + // ************** + // **** GRID **** + // ************** + + // create a linear grid with 20 cells + int cells = 20; + Grid grid = Grid(cells); + + MatrixXd concentrations = MatrixXd::Constant(1,20,20); + grid.setConcentrations(concentrations); + + + // ****************** + // **** BOUNDARY **** + // ****************** + + // create a boundary with constant values + Boundary bc = Boundary(grid, BC_TYPE_CONSTANT); + + + // ************************ + // **** SIMULATION ENV **** + // ************************ + + // set up a simulation environment + Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + + // (optional) set the timestep of the simulation + simulation.setTimestep(0.1); // timestep + + // (optional) set the number of iterations + simulation.setIterations(100); + + // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_OFF); + + // **** RUN SIMULATION **** + + // run the simulation + simulation.run(); +} \ No newline at end of file diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index 1490ac5..5f159bb 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -49,12 +49,12 @@ int main(int argc, char *argv[]) { // (optional) set boundary condition values for one side, e.g.: // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values - VectorXd bc_zero_values = VectorXd::Constant(20,0); - bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); - bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); - VectorXd bc_front_values = VectorXd::Constant(20,2000); - bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); - bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); + // VectorXd bc_zero_values = VectorXd::Constant(20,0); + // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); + // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); + // VectorXd bc_front_values = VectorXd::Constant(20,2000); + // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); + // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); // ************************ diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index d633503..df3ad0f 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -4,6 +4,7 @@ #include #include "Grid.hpp" +using namespace std; using namespace Eigen; enum BC_TYPE { @@ -18,6 +19,45 @@ enum BC_SIDE { BC_SIDE_BOTTOM }; +/******************************************* +class WallElement { + public: + WallElement() { + this->type = BC_TYPE_CLOSED; + this->value = 0; + } + + WallElement(double value) { + this->type = BC_TYPE_CONSTANT; + this->value = value; + } + + BC_TYPE getType() { + return this->type; + } + + double getValue() { + return this->value; + } + + private: + BC_TYPE type; + double value; +}; + +class BoundaryWall { + public: + BoundaryWall(int length) { + // create array with length many wall elements + } + + private: + BC_SIDE side; + int length; + vector wall; +}; +***********************/ + class Boundary { public: @@ -55,6 +95,10 @@ class Boundary { private: Grid grid; + + // need a way to save the bc type and value for each single 'boundary cell' + // perhaps an array for each side with structs containing the bc type as well as a value + // or another object that contains one boundary side BC_TYPE type; VectorXd left, right, top, bottom; }; diff --git a/proto/FTCS.ipynb b/proto/FTCS.ipynb index b64c934..4708428 100644 --- a/proto/FTCS.ipynb +++ b/proto/FTCS.ipynb @@ -458,7 +458,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.11.4" }, "orig_nbformat": 4 }, diff --git a/scripts/Adi2D_Reference.R b/scripts/Adi2D_Reference.R index 9b0bdd1..947ddc8 100644 --- a/scripts/Adi2D_Reference.R +++ b/scripts/Adi2D_Reference.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2023-01-05 17:52:55 delucia" +## Time-stamp: "Last modified 2023-05-11 17:31:41 delucia" ## Brutal implementation of 2D ADI scheme ## Square NxN grid with dx=dy=1 @@ -272,8 +272,8 @@ ADIHetDir <- function(field, dt, iter, alpha) { for (it in seq(1, iter)) { for (i in seq(2, ny-1)) { - Aij <- cbind(harm(alpha[i,], alpha[i-1,]), harm(alpha[i,], alpha[i+1,])) - Bij <- cbind(harm(alpha[,i], alpha[,i-1]), harm(alpha[,i], alpha[,i+1])) + Aij <- cbind(colMeans(rbind(alpha[i,], alpha[i-1,])), colMeans(rbind(alpha[i,], alpha[i+1,]))) + Bij <- cbind(rowMeans(cbind(alpha[,i], alpha[,i-1])), rowMeans(cbind(alpha[,i], alpha[,i+1]))) tmpX[i,] <- SweepByRowHetDir(i, res, dt=dt, Aij, Bij) } resY <- t(tmpX) @@ -321,22 +321,27 @@ SweepByRowHetDir <- function(i, field, dt, Aij, Bij) { ## adi2 <- ADI(n=51, dt=10, iter=200, alpha=1E-3) ## ref2 <- DoRef(n=51, alpha=1E-3, dt=10, iter=200) -n <- 5 +n <- 51 field <- matrix(0, n, n) alphas <- matrix(1E-5*runif(n*n, 1,2), n, n) -alphas1 <- matrix(3E-5, n, 25) -alphas2 <- matrix(1E-5, n, 26) -alphas <- cbind(alphas1, alphas2) +## dim(field) +## dim(alphas) +## all.equal(dim(field), dim(alphas)) + +## alphas1 <- matrix(3E-5, n, 25) +## alphas2 <- matrix(1E-5, n, 26) + +## alphas <- cbind(alphas1, alphas2) ## for (i in seq(1,nrow(alphas))) ## alphas[i,] <- seq(1E-7,1E-3, length=n) #diag(alphas) <- rep(1E-2, n) -adih <- ADIHetDir(field=field, dt=10, iter=200, alpha=alphas) -adi2 <- ADI(n=n, dt=10, iter=200, alpha=1E-5) +adih <- ADIHetDir(field=field, dt=20, iter=500, alpha=alphas) +adi2 <- ADI(n=n, dt=20, iter=500, alpha=1E-5) par(mfrow=c(1,3)) @@ -347,6 +352,17 @@ plot(adih[[length(adih)]], adi2[[length(adi2)]], pch=4, log="xy") abline(0,1) +cchet <- lapply(adih, round, digits=6) +cchom <- lapply(adi2, round, digits=6) + +plot(cchet[[length(cchet)]], cchom[[length(cchom)]], pch=4, log="xy", xlim=c(1e-6,1), ylim=c(1e-6,1)) +abline(0,1) + +cchet[[500]] + +str(adih) + + sapply(adih, sum) sapply(adi2, sum) @@ -360,3 +376,108 @@ image(adih[[length(adih)]]) points(0.5,0.5, col="red",pch=4) options(width=110) + + +FTCS_2D <- function(field, dt, iter, alpha) { + + if (!all.equal(dim(field), dim(alpha))) + stop("field and alpha are not matrix") + + ## now both field and alpha must be nx*ny matrices + nx <- ncol(field) + ny <- nrow(field) + dx <- dy <- 1 + + ## find out the center of the grid to apply conc=1 + cenx <- ceiling(nx/2) + ceny <- ceiling(ny/2) + field[cenx, ceny] <- 1 + + ## prepare containers for computations and outputs + tmp <- res <- field + + cflt <- 1/max(alpha)/4 + cat(":: CFL allowable time step: ", cflt,"\n") + + ## inner iterations + inner <- floor(dt/cflt) + if (inner == 0) { + ## dt < cflt, no inner iterations + inner <- 1 + tsteps <- dt + cat(":: No inner iter. required\n") + } else { + tsteps <- c(rep(cflt, inner), dt-inner*cflt) + cat(":: Number of inner iter. required: ", inner,"\n") + } + + + out <- vector(mode="list", length=iter) + + for (it in seq(1, iter)) { + cat(":: outer it: ", it) + + for (innerit in seq_len(inner)) { + for (i in seq(2, ny-1)) { + for (j in seq(2, nx-1)) { + ## tmp[i,j] <- res[i,j] + + ## + tsteps[innerit]/dx/dx * (res[i+1,j]*mean(alpha[i+1,j],alpha[i,j]) - + ## res[i,j] *(mean(alpha[i+1,j],alpha[i,j])+mean(alpha[i-1,j],alpha[i,j])) + + ## res[i-1,j]*mean(alpha[i-1,j],alpha[i,j])) + + ## + tsteps[innerit]/dy/dy * (res[i,j+1]*mean(alpha[i,j+1],alpha[i,j]) - + ## res[i,j] *(mean(alpha[i,j+1],alpha[i,j])+mean(alpha[i,j-1],alpha[i,j])) + + ## res[i,j-1]*mean(alpha[i,j-1],alpha[i,j])) + tmp[i,j] <- res[i,j] + + + tsteps[innerit]/dx/dx * ((res[i+1,j]-res[i,j]) * mean(alpha[i+1,j],alpha[i,j]) - + (res[i,j]-res[i-1,j]) * mean(alpha[i-1,j],alpha[i,j])) + + + tsteps[innerit]/dx/dx * ((res[i,j+1]-res[i,j]) * mean(alpha[i,j+1],alpha[i,j]) - + (res[i,j]-res[i,j-1]) * mean(alpha[i,j-1],alpha[i,j])) + } + } + ## swap back tmp to res for the next inner iteration + res <- tmp + } + cat("- done\n") + ## at end of inner it we store + out[[it]] <- res + } + + return(out) +} + +## testing that FTCS with homog alphas reverts to ADI/Reference sim +n <- 51 +field <- matrix(0, n, n) +alphas <- matrix(1E-3, n, n) + +adi2 <- ADI(n=51, dt=100, iter=20, alpha=1E-3) + +ref <- DoRef(n=51, alpha=1E-3, dt=100, iter=20) + +adihet <- ADIHetDir(field=field, dt=100, iter=20, alpha=alphas) + +ftcsh <- FTCS_2D(field=field, dt=100, iter=20, alpha=alphas) + + +par(mfrow=c(2,4)) +image(ref, main="Reference ODE.2D") +points(0.5,0.5, col="red",pch=4) +image(ftcsh[[length(ftcsh)]], main="FTCS 2D") +points(0.5,0.5, col="red",pch=4) +image(adihet[[length(adihet)]], main="ADI Heter.") +points(0.5,0.5, col="red",pch=4) +image(adi2[[length(adi2)]], main="ADI Homog.", col=terrain.colors(12)) +points(0.5,0.5, col="red",pch=4) +plot(ftcsh[[length(ftcsh)]], ref, pch=4, log="xy", xlim=c(1E-16, 1), ylim=c(1E-16, 1), + main = "FTCS_2D vs ref", xlab="FTCS 2D", ylab="Reference") +abline(0,1) +plot(ftcsh[[length(ftcsh)]], adihet[[length(adihet)]], pch=4, log="xy", xlim=c(1E-16, 1), ylim=c(1E-16, 1), + main = "FTCS_2D vs ADI Het", xlab="FTCS 2D", ylab="ADI 2D Heter.") +abline(0,1) +plot(ftcsh[[length(ftcsh)]], adi2[[length(adi2)]], pch=4, log="xy", xlim=c(1E-16, 1), ylim=c(1E-16, 1), + main = "FTCS_2D vs ADI Hom", xlab="FTCS 2D", ylab="ADI 2D Hom.") +abline(0,1) +plot(adihet[[length(adihet)]], adi2[[length(adi2)]], pch=4, log="xy", xlim=c(1E-16, 1), ylim=c(1E-16, 1), + main = "ADI Het vs ADI Hom", xlab="ADI Het", ylab="ADI 2D Hom.") +abline(0,1) + diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 6db41de..58812a9 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -4,7 +4,7 @@ using namespace std; -double calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = false) { +double calcAlphaIntercell(double alpha1, double alpha2, bool useHarmonic = false) { if (useHarmonic) { return 2 / ((1/alpha1) + (1/alpha2)); } else { @@ -12,7 +12,168 @@ double calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = fal } } -MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) { + +double calcHorizontalChange(Grid grid, int row, int col) { + + double result = + calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) + * grid.getConcentrations()(row,col+1) + - ( + calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) + + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) + ) + * grid.getConcentrations()(row,col) + + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) + * grid.getConcentrations()(row,col-1); + + return result; +} + + +double calcVerticalChange(Grid grid, int row, int col) { + + double result = + calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) + * grid.getConcentrations()(row+1,col) + - ( + calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) + + calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)) + ) + * grid.getConcentrations()(row,col) + + calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)) + * grid.getConcentrations()(row-1,col); + + return result; +} + + +double calcHorizontalChangeLeftBoundaryConstant(Grid grid, Boundary bc, int row, int col) { + + double result = + calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) + * grid.getConcentrations()(row,col+1) + - ( + calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) + + 2 * grid.getAlphaX()(row,col) + ) + * grid.getConcentrations()(row,col) + + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row); + + return result; +} + + +double calcHorizontalChangeLeftBoundaryClosed() { + return 0; +} + + +double calcHorizontalChangeRightBoundaryConstant(Grid grid, Boundary bc, int row, int col) { + + double result = + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row) + - ( + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) + + 2 * grid.getAlphaX()(row,col) + ) + * grid.getConcentrations()(row,col) + + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) + * grid.getConcentrations()(row,col-1); + + return result; +} + + +double calcHorizontalChangeRightBoundaryClosed() { + return 0; +} + + +double calcVerticalChangeTopBoundaryConstant(Grid grid, Boundary bc, int row, int col) { + + double result = + calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col)) + * grid.getConcentrations()(row+1,col) + - ( + calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col)) + + 2 * grid.getAlphaY()(row, col) + ) + * grid.getConcentrations()(row, col) + + 2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col); + + return result; +} + + +double calcVerticalChangeTopBoundaryClosed() { + return 0; +} + + +double calcVerticalChangeBottomBoundaryConstant(Grid grid, Boundary bc, int row, int col) { + + double result = + 2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col) + - ( + calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) + + 2 * grid.getAlphaY()(row, col) + ) + * grid.getConcentrations()(row, col) + + calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) + * grid.getConcentrations()(row-1,col); + + return result; +} + + +double calcVerticalChangeBottomBoundaryClosed() { + return 0; +} + + +MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) { + int colMax = grid.getCol(); + double deltaCol = grid.getDeltaCol(); + + MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0); + + // only one row in 1D case + int row = 0; + + // inner cells + for (int col = 1; col < colMax-1; col++) { + concentrations_t1(row,col) = grid.getConcentrations()(row,col) + + timestep / (deltaCol*deltaCol) + * ( + calcHorizontalChange(grid, row, col) + ) + ; + } + + // left boundary + int col = 0; + concentrations_t1(row, col) = grid.getConcentrations()(row,col) + + timestep / (deltaCol*deltaCol) + * ( + calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col) + ) + ; + + + // right boundary + col = colMax-1; + concentrations_t1(row,col) = grid.getConcentrations()(row,col) + + timestep / (deltaCol*deltaCol) + * ( + calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col) + ) + ; + + return concentrations_t1; +} + + +MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) { int rowMax = grid.getRow(); int colMax = grid.getCol(); double deltaRow = grid.getDeltaRow(); @@ -20,32 +181,20 @@ MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) { // Matrix with concentrations at time t+1 // TODO profiler / only use 2 matrices - MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 1); + MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); // inner cells - // (should have 7 calls to current concentration) + // these do not depend on the boundary condition type for (int row = 1; row < rowMax-1; row++) { for (int col = 1; col < colMax-1; col++) { concentrations_t1(row, col) = grid.getConcentrations()(row, col) + timestep / (deltaRow*deltaRow) * ( - calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) - * grid.getConcentrations()(row+1,col) - - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) - + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))) - * grid.getConcentrations()(row,col) - + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)) - * grid.getConcentrations()(row-1,col) + calcVerticalChange(grid, row, col) ) + timestep / (deltaCol*deltaCol) * ( - calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col+1) - - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))) - * grid.getConcentrations()(row,col) - + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col-1) + calcHorizontalChange(grid, row, col) ) ; } @@ -53,184 +202,119 @@ MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) { // boundary conditions // left without corners / looping over rows - // (should have 6 calls to current concentration) int col = 0; for (int row = 1; row < rowMax-1; row++) { concentrations_t1(row, col) = grid.getConcentrations()(row,col) + timestep / (deltaCol*deltaCol) - * (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col+1) - - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - + 2 * grid.getAlphaX()(row,col)) * grid.getConcentrations()(row,col) - + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row)) + * ( + calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col) + ) + timestep / (deltaRow*deltaRow) - * (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) - * grid.getConcentrations()(row+1,col) - - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) - + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))) - * grid.getConcentrations()(row,col) - + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)) - * grid.getConcentrations()(row-1,col)); + * ( + calcVerticalChange(grid, row, col) + ) + ; } // right without corners / looping over columns - // (should have 6 calls to current concentration) col = colMax-1; for (int row = 1; row < rowMax-1; row++) { concentrations_t1(row,col) = grid.getConcentrations()(row,col) + timestep / (deltaCol*deltaCol) - * (2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row) - - (calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) - + 2 * grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col) - + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col-1)) + * ( + calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col) + ) + timestep / (deltaRow*deltaRow) - * (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) - * grid.getConcentrations()(row+1,col) - - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) - + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))) - * grid.getConcentrations()(row,col) - + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)) - * grid.getConcentrations()(row-1,col)); + * ( + calcVerticalChange(grid, row, col) + ) + ; } // top without corners / looping over cols - // (should have 6 calls to current concentration) int row = 0; for (int col=1; col this check is last in order + // - } diff --git a/src/Simulation.cpp b/src/Simulation.cpp index a078f47..6d7bc6f 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -84,7 +84,11 @@ void Simulation::run() { } printConcentrationsConsole(); if (csv_output >= CSV_OUTPUT_ON) { - printConcentrationsCSV("test", true); + bool append = false; + if (csv_output == CSV_OUTPUT_VERBOSE) { + append = true; + } + printConcentrationsCSV("test", append); } } else if (approach == BTCS_APPROACH) { for (int i = 0; i < iterations; i++) { diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1096dd3..ed46443 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,17 +1,20 @@ -include(FetchContent) +find_library(DOCTEST_LIB doctest) -FetchContent_Declare( +if(NOT DOCTEST_LIB) + include(FetchContent) + + FetchContent_Declare( DocTest GIT_REPOSITORY https://github.com/doctest/doctest.git - GIT_TAG v2.4.9 -) + GIT_TAG v2.4.9) -FetchContent_MakeAvailable(DocTest) + FetchContent_MakeAvailable(DocTest) +endif() add_executable(testTug setup.cpp testBoundaryCondition.cpp testDiffusion.cpp) target_link_libraries(testTug doctest tug) -add_custom_target(check - COMMAND $ - DEPENDS testTug -) +add_custom_target( + check + COMMAND $ + DEPENDS testTug) diff --git a/utils/ci.Dockerfile b/utils/ci.Dockerfile new file mode 100644 index 0000000..30ddc7a --- /dev/null +++ b/utils/ci.Dockerfile @@ -0,0 +1,4 @@ +FROM alpine + MAINTAINER Max Luebke + +RUN apk add --no-cache build-base openmp cmake git eigen-dev