From 32b05a8a871b735eb817790cf0015384818d32b9 Mon Sep 17 00:00:00 2001 From: philippun Date: Wed, 23 Aug 2023 12:24:35 +0200 Subject: [PATCH] added Thomas Solver with option to choose solver and cleaned up the repository --- examples/CMakeLists.txt | 15 +- examples/Rcpp-BTCS-1d.cpp | 59 ------ examples/Rcpp-BTCS-2d.cpp | 60 ------ examples/Rcpp-interface.R | 159 --------------- examples/Rcpp-interface.cpp | 59 ------ examples/RcppFTCS.cpp | 40 ---- examples/boundary_example1D.cpp | 35 ---- examples/first_example.cpp | 72 ------- examples/main_1D.cpp | 55 ----- examples/main_2D.cpp | 58 ------ examples/main_2D_const.cpp | 59 ------ examples/main_2D_mdl.cpp | 67 ------ examples/meeting_example.cpp | 38 ---- examples/profiling.cpp | 50 ----- examples/profiling_openmp.cpp | 5 +- examples/pseudo_example.cpp | 25 --- examples/second_example.cpp | 93 --------- include/tug/BoundaryCondition.hpp | 283 -------------------------- include/tug/Diffusion.hpp | 140 ------------- include/tug/Simulation.hpp | 21 +- include/tug/Solver.hpp | 41 ---- include/tug/progressbar.hpp | 191 ----------------- src/BTCS.cpp | 328 ------------------------------ src/BTCSv2.cpp | 76 ++++++- src/Boundary.cpp | 3 +- src/BoundaryCondition.cpp | 209 ------------------- src/CMakeLists.txt | 2 +- src/README.md | 9 - src/Simulation.cpp | 132 +++++++----- src/Solver.cpp | 61 ------ test/CMakeLists.txt | 2 +- test/testBoundaryCondition.cpp | 196 ------------------ test/testDiffusion.cpp | 170 ---------------- 33 files changed, 181 insertions(+), 2632 deletions(-) delete mode 100644 examples/Rcpp-BTCS-1d.cpp delete mode 100644 examples/Rcpp-BTCS-2d.cpp delete mode 100644 examples/Rcpp-interface.R delete mode 100644 examples/Rcpp-interface.cpp delete mode 100644 examples/RcppFTCS.cpp delete mode 100644 examples/boundary_example1D.cpp delete mode 100644 examples/first_example.cpp delete mode 100644 examples/main_1D.cpp delete mode 100644 examples/main_2D.cpp delete mode 100644 examples/main_2D_const.cpp delete mode 100644 examples/main_2D_mdl.cpp delete mode 100644 examples/meeting_example.cpp delete mode 100644 examples/profiling.cpp delete mode 100644 examples/pseudo_example.cpp delete mode 100644 examples/second_example.cpp delete mode 100644 include/tug/BoundaryCondition.hpp delete mode 100644 include/tug/Diffusion.hpp delete mode 100644 include/tug/Solver.hpp delete mode 100644 include/tug/progressbar.hpp delete mode 100644 src/BTCS.cpp delete mode 100644 src/BoundaryCondition.cpp delete mode 100644 src/Solver.cpp delete mode 100644 test/testBoundaryCondition.cpp delete mode 100644 test/testDiffusion.cpp diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 33aeea2..91c998a 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,24 +1,17 @@ -add_executable(first_example first_example.cpp) -add_executable(second_example second_example.cpp) -add_executable(boundary_example1D boundary_example1D.cpp) add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp) add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp) add_executable(BTCS_1D_proto_example BTCS_1D_proto_example.cpp) add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp) -add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp) add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp) -add_executable(meeting_example meeting_example.cpp) -target_link_libraries(first_example tug) -target_link_libraries(second_example tug) -target_link_libraries(boundary_example1D tug) +add_executable(profiling_openmp profiling_openmp.cpp) target_link_libraries(FTCS_1D_proto_example tug) target_link_libraries(FTCS_2D_proto_example tug) target_link_libraries(BTCS_1D_proto_example tug) target_link_libraries(BTCS_2D_proto_example tug) -target_link_libraries(FTCS_2D_proto_example_mdl tug) target_link_libraries(reference-FTCS_2D_closed tug) -target_link_libraries(meeting_example tug) -# target_link_libraries(FTCS_2D_proto_example easy_profiler) +target_link_libraries(profiling_openmp tug) +add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp) add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp) target_link_libraries(FTCS_2D_proto_closed_mdl tug) +target_link_libraries(FTCS_2D_proto_example_mdl tug) diff --git a/examples/Rcpp-BTCS-1d.cpp b/examples/Rcpp-BTCS-1d.cpp deleted file mode 100644 index 8949d85..0000000 --- a/examples/Rcpp-BTCS-1d.cpp +++ /dev/null @@ -1,59 +0,0 @@ -#include "../include/diffusion/BTCSDiffusion.hpp" -#include "../include/diffusion/BoundaryCondition.hpp" -#include // for copy, max -#include -#include -#include // for std -#include // for vector -#include - -using namespace std; -using namespace Diffusion; -using namespace Rcpp; -//using namespace Eigen; - -// [[Rcpp::depends(RcppEigen)]] -// [[Rcpp::plugins("cpp11")]] -// [[Rcpp::export]] -std::vector & diff1D(int n, - double length, - std::vector & field, - std::vector & alpha, - double timestep, - double bc_left, - double bc_right, - int iterations) { - // dimension of grid - int dim = 1; - - // 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); - - diffu.setXDimensions(length, n); - - // set the boundary condition for the left ghost cell to dirichlet - bc[0] = {Diffusion::BC_CONSTANT, bc_left}; - bc[n-1] = {Diffusion::BC_CONSTANT, bc_right}; - - // set timestep for simulation to 1 second - diffu.setTimestep(timestep); - - //cout << setprecision(12); - - // loop 100 times - // output is currently generated by the method itself - for (int i = 0; i < iterations; i++) { - diffu.simulate(field.data(), alpha.data(), bc.data()); - } - - // for (auto & cell : field) { - // Rcout << cell << "\n"; - // } - - return(field); -} diff --git a/examples/Rcpp-BTCS-2d.cpp b/examples/Rcpp-BTCS-2d.cpp deleted file mode 100644 index 3a5e404..0000000 --- a/examples/Rcpp-BTCS-2d.cpp +++ /dev/null @@ -1,60 +0,0 @@ -#include "../include/diffusion/BTCSDiffusion.hpp" -#include "../include/diffusion/BoundaryCondition.hpp" -#include // for copy, max -#include -#include -#include // for std -#include // for vector -#include - -using namespace std; -using namespace Diffusion; -using namespace Rcpp; -//using namespace Eigen; - -// [[Rcpp::depends(RcppEigen)]] -// [[Rcpp::plugins("cpp11")]] -// [[Rcpp::export]] -std::vector & diff2D(int nx, - int ny, - double lenx, - double leny, - std::vector & field, - std::vector & alpha, - double timestep, - int iterations) -{ - // problem dimensionality - int dim = 2; - // total number of grid cells - int n = nx*ny; - - std::vector bc(n, {0,0}); - - // create instance of diffusion module - BTCSDiffusion diffu(dim); - - diffu.setXDimensions(lenx, nx); - diffu.setXDimensions(leny, ny); - - // set the boundary condition for the left ghost cell to dirichlet - // bc[0] = {Diffusion::BC_CONSTANT, bc_left}; - // bc[n-1] = {Diffusion::BC_CONSTANT, bc_right}; - - // set timestep for simulation to 1 second - diffu.setTimestep(timestep); - - //cout << setprecision(12); - - // loop 100 times - // output is currently generated by the method itself - for (int i = 0; i < iterations; i++) { - diffu.simulate(field.data(), alpha.data(), bc.data()); - } - - // for (auto & cell : field) { - // Rcout << cell << "\n"; - // } - - return(field); -} diff --git a/examples/Rcpp-interface.R b/examples/Rcpp-interface.R deleted file mode 100644 index fc1280b..0000000 --- a/examples/Rcpp-interface.R +++ /dev/null @@ -1,159 +0,0 @@ -## Time-stamp: "Last modified 2022-03-16 14:01:11 delucia" -library(Rcpp) -library(RcppEigen) -library(ReacTran) -library(deSolve) - -options(width=110) - -setwd("app") - -## This creates the "diff1D" function with our BTCSdiffusion code -sourceCpp("Rcpp-BTCS-1d.cpp") - -### FTCS explicit (same name) -sourceCpp("RcppFTCS.cpp") - - - -## Grid 101 -## Set initial conditions -N <- 1001 -D.coeff <- 1E-3 -C0 <- 1 ## Initial concentration (mg/L) -X0 <- 0 ## Location of initial concentration (m) -## Yini <- c(C0, rep(0,N-1)) - -## Ode1d solution -xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = N) -x <- xgrid$x.mid -Diffusion <- function (t, Y, parms){ - tran <- tran.1D(C = Y, C.up = 0, C.down = 0, D = parms$D, dx = xgrid) - return(list(tran$dC)) -} - - -## gaussian pulse as initial condition -sigma <- 0.02 -Yini <- 0.5*exp(-0.5*((x-1/2.0)**2)/sigma**2) - -## plot(x, Yini, type="l") - -parms1 <- list(D=D.coeff) -# 1 timestep, 10 s -times <- seq(from = 0, to = 1, by = 0.1) - -system.time({ - out1 <- ode.1D(y = Yini, times = times, func = Diffusion, - parms = parms1, dimens = N)[11,-1] -}) - -## Now with BTCS -alpha <- rep(D.coeff, N) - -system.time({ - out2 <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = 0.1, 0, 0, iterations = 10) -}) - -## plot(out1, out2) -## abline(0,1) - -## matplot(cbind(out1,out2),type="l", col=c("black","red"),lty="solid", lwd=2, -## xlab="grid element", ylab="Concentration", las=1) -## legend("topright", c("ReacTran ode1D", "BTCS 1d"), text.col=c("black","red"), bty = "n") - - - - - -system.time({ - out3 <- RcppFTCS(n=N, length=1, field=Yini, alpha=1E-3, bc_left = 0, bc_right = 0, timestep = 1) -}) - -## Poor man's -mm <- colMeans(rbind(out2,out3)) - - - -matplot(cbind(Yini,out1, out2, out3, mm),type="l", col=c("grey","black","red","blue","green4"), lty="solid", lwd=2, - xlab="grid element", ylab="Concentration", las=1) -legend("topright", c("init","ReacTran ode1D", "BTCS 1d", "FTCS", "poor man's CN"), text.col=c("grey","black","red","blue","green4"), bty = "n") - - -sum(Yini) -sum(out1) -sum(out2) -sum(out3) -sum(mm) - -## Yini <- 0.2*sin(pi/0.1*x)+0.2 -## plot(Yini) - -## plot(out3) - -Fun <- function(dx) { - tmp <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = floor(1/dx)) - sqrt(sum({out1-tmp}^2)) -} - -reso <- optimise(f=Fun, interval=c(1E-5, 1E-1), maximum = FALSE) - - -dx <- 0.0006038284 -floor(1/dx) - -1/dx - -system.time({ - out2o <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = 1656) -}) - -matplot(cbind(out1, out2o),type="l", col=c("black","red"), lty="solid", lwd=2, - xlab="grid element", ylab="Concentration", las=1) -legend("topright", c("ReacTran ode1D", "BTCS 1d dx=0.0006"), text.col=c("black","red"), bty = "n") - - -dx <- 0.05 - -system.time({ - out2o <- diff1D(n=N, length=1, field=Yini, alpha=alpha, timestep = dx, 0, 0, iterations = 1/dx) -}) - -matplot(cbind(out1, out2o),type="l", col=c("black","red"), lty="solid", lwd=2, - xlab="grid element", ylab="Concentration", las=1) -legend("topright", c("ReacTran ode1D", "BTCS 1d dx=0.0006"), text.col=c("black","red"), bty = "n") - -Matplot - - - - -## This creates the "diff1D" function with our BTCSdiffusion code -sourceCpp("Rcpp-BTCS-2d.cpp") - -n <- 256 -a2d <- rep(1E-3, n^2) - -init2d <- readRDS("gs1.rds") - -ll <- {init2d - min(init2d)}/diff(range(init2d)) - -system.time({ - res1 <- diff2D(nx=N, ny=N, lenx=1, leny=1, field=ll, alpha=a2d, timestep = 0.1, iterations = 10) -}) - -hist(ll,32) - - - - - - - - - - - - - - diff --git a/examples/Rcpp-interface.cpp b/examples/Rcpp-interface.cpp deleted file mode 100644 index 8949d85..0000000 --- a/examples/Rcpp-interface.cpp +++ /dev/null @@ -1,59 +0,0 @@ -#include "../include/diffusion/BTCSDiffusion.hpp" -#include "../include/diffusion/BoundaryCondition.hpp" -#include // for copy, max -#include -#include -#include // for std -#include // for vector -#include - -using namespace std; -using namespace Diffusion; -using namespace Rcpp; -//using namespace Eigen; - -// [[Rcpp::depends(RcppEigen)]] -// [[Rcpp::plugins("cpp11")]] -// [[Rcpp::export]] -std::vector & diff1D(int n, - double length, - std::vector & field, - std::vector & alpha, - double timestep, - double bc_left, - double bc_right, - int iterations) { - // dimension of grid - int dim = 1; - - // 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); - - diffu.setXDimensions(length, n); - - // set the boundary condition for the left ghost cell to dirichlet - bc[0] = {Diffusion::BC_CONSTANT, bc_left}; - bc[n-1] = {Diffusion::BC_CONSTANT, bc_right}; - - // set timestep for simulation to 1 second - diffu.setTimestep(timestep); - - //cout << setprecision(12); - - // loop 100 times - // output is currently generated by the method itself - for (int i = 0; i < iterations; i++) { - diffu.simulate(field.data(), alpha.data(), bc.data()); - } - - // for (auto & cell : field) { - // Rcout << cell << "\n"; - // } - - return(field); -} diff --git a/examples/RcppFTCS.cpp b/examples/RcppFTCS.cpp deleted file mode 100644 index dbe473b..0000000 --- a/examples/RcppFTCS.cpp +++ /dev/null @@ -1,40 +0,0 @@ -// Time-stamp: "Last modified 2022-03-15 18:15:39 delucia" -#include -#include // for std -#include // for vector - -using namespace std; -using namespace Rcpp; - -// [[Rcpp::plugins("cpp11")]] -// [[Rcpp::export]] -NumericVector RcppFTCS(int n, - double length, - NumericVector & field, - double alpha, - double bc_left, - double bc_right, - double timestep) -{ - // dimension of grid - NumericVector ext (clone(field)); - double dx = length / ((double) n - 1.); - double dt = 0.25*dx*dx/alpha; - - - double afac = alpha*dt/dx/dx; - int iter = (int) (timestep/dt); - - Rcout << "dt: " << dt << "; inner iterations: " << iter << endl; - - - for (int it = 0; it < iter; it++){ - for (int i = 1; i < ext.size()-1; i++) { - ext[i] = (1. - 2*afac)*ext[i] + afac*(ext[i+1]+ext[i-1]); - } - ext[0] = bc_left; - ext[n-1] = bc_right; - } - return(ext); -} - diff --git a/examples/boundary_example1D.cpp b/examples/boundary_example1D.cpp deleted file mode 100644 index 4237f11..0000000 --- a/examples/boundary_example1D.cpp +++ /dev/null @@ -1,35 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace tug::bc; -using namespace tug::diffusion; - -int main(int argc, char *argv[]) { - TugInput input; - - BoundaryCondition example(10); - - uint8_t side = BC_SIDE_LEFT; - boundary_condition bc; - bc.type = BC_TYPE_CONSTANT; - bc.value = 1; - input.setBoundaryCondition(example); - BoundaryCondition returnedBC = input.getBoundaryCondition(); - - // example.setSide(side, bc); - vector result_left = example.getSide(side); - // vector result_top = example.getSide(BC_SIDE_TOP); - - for (auto i : result_left) { - cout << i.value << " "; - } - cout << endl; - // for (auto i : result_top) { - // cout << i.value << " "; - // } - // cout << endl; -} \ No newline at end of file diff --git a/examples/first_example.cpp b/examples/first_example.cpp deleted file mode 100644 index e3c3549..0000000 --- a/examples/first_example.cpp +++ /dev/null @@ -1,72 +0,0 @@ -#include "tug/BoundaryCondition.hpp" -#include - -#include -#include -#include - -using namespace std; -using namespace tug::diffusion; -using namespace tug::bc; - -int main(int argc, char *argv[]) { - - // define problem dimensionality - // set grid sizes for each dimension - int dim = 1; - int n = 20; - - vector alpha(n, 1e-1); - // double alpha = 1e-1; - vector field(n, 1e-6); - for (int i = 1; i<20; i++) { - field[i] = 0; - } - // double field = 1e-6; - - - TugGrid grid_param; // why is grid_param defined separately? - // grid_param.grid_cells[0] = 20; - // grid_param.grid_cells[1] = 0; - // grid_param.grid_cells[2] = 0; - - // grid_param.domain_size[0] = 20; - // grid_param.domain_size[1] = 0; - // grid_param.domain_size[2] = 0; - - - - TugInput input_param; - input_param.setTimestep(1.); - //input_param.grid = grid_param; - input_param.setGridCellN(n); - input_param.setDomainSize(n); // what is domain????? - BoundaryCondition bc(n); - input_param.setBoundaryCondition(bc); - - BoundaryCondition bc2 = input_param.getBoundaryCondition(); - auto [bc_left, bc_right] = bc2.row_boundary(0); - cout << "left: " << unsigned(bc_left.type) << endl; - cout << "right: " << unsigned(bc_right.type) << endl; - - ofstream myfile; - myfile.open("output.csv"); - if (!myfile) { - exit(1); - } - - for (int t = 0; t < 10000; t++) { - double result = BTCS_1D(input_param, &field[0], &alpha[0]); - //myfile << result; - //myfile << '\n'; - //myfile << "Vector contents: "; - for (int i = 0; i < field.size(); i++) { - myfile << field[i]; - if (i < field.size()-1) { - myfile << ", "; - } - } - myfile << std::endl; - } - myfile.close(); -} \ No newline at end of file diff --git a/examples/main_1D.cpp b/examples/main_1D.cpp deleted file mode 100644 index 1b143fc..0000000 --- a/examples/main_1D.cpp +++ /dev/null @@ -1,55 +0,0 @@ -#include "diffusion/BTCSBoundaryCondition.hpp" -#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 = 1; - - int n = 20; - - // create input + diffusion coefficients for each grid cell - std::vector alpha(n, 1e-1); - std::vector field(n, 1e-6); - - BTCSBoundaryCondition bc(n); - - // create instance of diffusion module - BTCSDiffusion diffu(dim); - - diffu.setXDimensions(1, n); - - // set the boundary condition for the left ghost cell to dirichlet - bc(BC_SIDE_LEFT) = {BC_TYPE_CONSTANT, 5e-6}; - // bc[0] = {Diffusion::BC_CONSTANT, 5e-6}; - // diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT, - // 5. * std::pow(10, -6)); - - // set timestep for simulation to 1 second - diffu.setTimestep(1.); - - cout << setprecision(12); - - // loop 100 times - // output is currently generated by the method itself - for (int i = 0; i < 100; i++) { - diffu.simulate(field.data(), alpha.data(), bc); - - cout << "Iteration: " << i << "\n\n"; - - for (auto &cell : field) { - cout << cell << "\n"; - } - - cout << "\n" << endl; - } - - return 0; -} diff --git a/examples/main_2D.cpp b/examples/main_2D.cpp deleted file mode 100644 index c76737d..0000000 --- a/examples/main_2D.cpp +++ /dev/null @@ -1,58 +0,0 @@ -#include "diffusion/BTCSBoundaryCondition.hpp" -#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, 1e-6); - std::vector field(n * m, 0); - BTCSBoundaryCondition bc(n, m); - - // create instance of diffusion module - BTCSDiffusion diffu(dim); - - diffu.setXDimensions(n, n); - diffu.setYDimensions(m, m); - - // set inlet to higher concentration without setting bc - field[12] = 1; - - // set timestep for simulation to 1 second - diffu.setTimestep(1); - - cout << setprecision(12); - - for (int t = 0; t < 1000; t++) { - diffu.simulate(field.data(), alpha.data(), bc); - - cout << "Iteration: " << t << "\n\n"; - - double sum = 0; - - // 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]; - sum += field[i * n + j]; - } - cout << "\n"; - } - - cout << "sum: " << sum << "\n" << endl; - } - - return 0; -} diff --git a/examples/main_2D_const.cpp b/examples/main_2D_const.cpp deleted file mode 100644 index a9a2a29..0000000 --- a/examples/main_2D_const.cpp +++ /dev/null @@ -1,59 +0,0 @@ -#include "diffusion/BTCSBoundaryCondition.hpp" -#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, 1e-1); - std::vector field(n * m, 1e-6); - BTCSBoundaryCondition bc(n, m); - - // create instance of diffusion module - BTCSDiffusion diffu(dim); - - diffu.setXDimensions(1, n); - diffu.setYDimensions(1, m); - - boundary_condition input = {Diffusion::BC_TYPE_CONSTANT, 5e-6}; - - bc.setSide(BC_SIDE_LEFT, input); - - // for (int i = 1; i <= n; i++) { - // bc[(n + 2) * i] = {Diffusion::BC_CONSTANT, 5e-6}; - // } - // 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); - - 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/examples/main_2D_mdl.cpp b/examples/main_2D_mdl.cpp deleted file mode 100644 index 16805af..0000000 --- a/examples/main_2D_mdl.cpp +++ /dev/null @@ -1,67 +0,0 @@ -#include "diffusion/BTCSBoundaryCondition.hpp" -#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, 1e-3); - std::vector field(n * m, 0.); - BTCSBoundaryCondition bc(n, m); - - 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++) { - double time = diffu.simulate(field.data(), alpha.data(), bc); - - cerr << "time elapsed: " << time << " seconds" << endl; - - cout << t; - - for (int i = 0; i < m * n; i++) { - cout << "," << field[i]; - } - cout << endl; - } - - return 0; -} diff --git a/examples/meeting_example.cpp b/examples/meeting_example.cpp deleted file mode 100644 index 27d5540..0000000 --- a/examples/meeting_example.cpp +++ /dev/null @@ -1,38 +0,0 @@ -#include "tug/Boundary.hpp" -#include - -int main(int argc, char *argv[]) { - int row = 20; - int col = 30; - - // grid - Grid grid = Grid(row,col); - MatrixXd concentrations = MatrixXd::Constant(row,col,0); - concentrations(10,4) = 2000; - grid.setConcentrations(concentrations); - MatrixXd alphaX = MatrixXd::Constant(row,col,1); - for (int i = 0; i < row; i++) { - alphaX(i,0) = 0.05; - } - MatrixXd alphaY = MatrixXd::Constant(row,col,1); - grid.setAlpha(alphaX, alphaY); - - // boundary - Boundary bc = Boundary(grid); - bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); - bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); - bc.setBoundarySideConstant(BC_SIDE_TOP, 0); - bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); - bc.setBoundaryElementClosed(BC_SIDE_LEFT, 3); - bc.setBoundaryElementClosed(BC_SIDE_LEFT, 4); - bc.setBoundaryElementClosed(BC_SIDE_LEFT, 5); - bc.setBoundaryElementConstant(BC_SIDE_BOTTOM, 17, 100); - - // simulation - Simulation sim = Simulation(grid, bc, BTCS_APPROACH); - sim.setTimestep(0.1); - sim.setIterations(100); - sim.setOutputCSV(CSV_OUTPUT_XTREME); - - sim.run(); -} \ No newline at end of file diff --git a/examples/profiling.cpp b/examples/profiling.cpp deleted file mode 100644 index a162705..0000000 --- a/examples/profiling.cpp +++ /dev/null @@ -1,50 +0,0 @@ -#include -#include -#include - -int main(int argc, char *argv[]) { - - int n[4] = {10, 20, 50, 100}; - int iterations[1] = {100}; - int repetition = 10; - - ofstream myfile; - myfile.open("btcs_time_measurement_openmp_10.csv"); - - for (int i = 0; i < size(n); i++){ - cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl; - myfile << "Grid size: " << n[i] << " x " << n[i] << endl; - for(int j = 0; j < size(iterations); j++){ - cout << "Iterations: " << iterations[j] << endl; - for (int k = 0; k < repetition; k++){ - cout << "Wiederholung: " << k << endl; - Grid grid = Grid(n[i], n[i]); - grid.setDomain(n[i], n[i]); - - MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); - concentrations(5,5) = 1; - grid.setConcentrations(concentrations); - MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5); - - Boundary bc = Boundary(grid); - - Simulation sim = Simulation(grid, bc, BTCS_APPROACH); - - - sim.setTimestep(0.001); - sim.setIterations(iterations[j]); - sim.setOutputCSV(CSV_OUTPUT_OFF); - - auto begin = std::chrono::high_resolution_clock::now(); - sim.run(); - auto end = std::chrono::high_resolution_clock::now(); - auto milliseconds = std::chrono::duration_cast(end - begin); - myfile << milliseconds.count() << endl; - } - } - cout << endl; - myfile << endl; - - } - myfile.close(); -} \ No newline at end of file diff --git a/examples/profiling_openmp.cpp b/examples/profiling_openmp.cpp index 9e62472..f41a0a0 100644 --- a/examples/profiling_openmp.cpp +++ b/examples/profiling_openmp.cpp @@ -9,7 +9,7 @@ int main(int argc, char *argv[]) { int repetition = 10; ofstream myfile; - myfile.open("time_measure_experiment_openmp_thread_6.csv"); + myfile.open("time_measure_experiment_openmp_thread_1_EigenLU.csv"); for (int i = 0; i < size(n); i++){ cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl; @@ -29,7 +29,8 @@ int main(int argc, char *argv[]) { Boundary bc = Boundary(grid); - Simulation sim = Simulation(grid, bc, FTCS_APPROACH); + Simulation sim = Simulation(grid, bc, BTCS_APPROACH); + sim.setSolver(EIGEN_LU_SOLVER); sim.setTimestep(0.001); diff --git a/examples/pseudo_example.cpp b/examples/pseudo_example.cpp deleted file mode 100644 index 09dfef2..0000000 --- a/examples/pseudo_example.cpp +++ /dev/null @@ -1,25 +0,0 @@ - -// create a new grid -// dimension is implicitly determined by number of arguments -// grid = new Grid(m=10, n=10) - -// array_x = array[] -// array_y = array[] -// grid.setAlpha_x(&array_x) -// grid.setAlpha_y(&array_y) - -// float time -// grid.setTimestep(time) - -// bc = new BoundaryCondition() -// grid.setBoundaryCondition(bc) - - - -// create new simulation run -// int iterations -// bool csv = true -// bool explicit = true -// simulation = new Simulation(grid, iterations) -// simulation.setSolver(BTCS, Thomas) -// simulation.run() \ No newline at end of file diff --git a/examples/second_example.cpp b/examples/second_example.cpp deleted file mode 100644 index e146cba..0000000 --- a/examples/second_example.cpp +++ /dev/null @@ -1,93 +0,0 @@ -#include "tug/BoundaryCondition.hpp" -#include - -#include -#include -#include - -using namespace std; -using namespace tug::diffusion; -using namespace tug::bc; - -int main(int argc, char *argv[]) { - - int dim = 2; - int n = 20; - int m = 20; - - vector alpha(n * m, 1); - vector field(n * m, 0); - field[0] = 2000; - // field[n * 19] = 2000; - // field[n * 19 + 19] = 2000; - - // for (int i = 1; i<20; i++) { - // for (int j = 0; j<20; j++ ) { - // field[i] = 0; - // } - // } - - // print field - cout << "Initial field:" << endl; - for (int i = 0; i -#include -#include -#include -#include - -typedef uint8_t bctype; - -namespace tug { -namespace bc { - -enum { - // flux = Neumann? Cauchy combination of Neumann and Dirichlet? - BC_TYPE_CLOSED, /**< Defines a closed boundary condition (special case of Neumann condition). */ - BC_TYPE_FLUX, /**< Defines a flux/Cauchy boundary condition. */ - BC_TYPE_CONSTANT, /**< Defines a constant/Dirichlet boundary condition. */ - BC_UNSET /**< Indicates undefined boundary condition*/ -}; - -enum { - BC_SIDE_LEFT, /**< Defines boundary conditions for the left side of the grid. - */ - BC_SIDE_RIGHT, /**< Defines boundary conditions for the right side of the - grid. */ - BC_SIDE_TOP, /**< Defines boundary conditions for the top of the grid. */ - BC_SIDE_BOTTOM, /**< Defines boundary conditions for the bottom of the grid. - */ - BC_INNER -}; - -/** - * Defines the boundary condition type and according value. - * QUESTION: For what is the struct bc necessary? - */ -typedef struct boundary_condition_s { - bctype type; /**< Type of the boundary condition */ - double value; /**< Value of the boundary condition. Either a concrete value of - concentration for BC_TYPE_CONSTANT or gradient when type is - BC_TYPE_FLUX. Unused if BC_TYPE_CLOSED.*/ -} boundary_condition; - -/** - * Represents both boundary conditions of a row/column. - */ -typedef std::array bc_tuple; -typedef std::vector bc_vec; - -/** - * Class to define the boundary condition of a grid. - */ -class BoundaryCondition { -public: - /** - * Creates a new instance with two elements. Used when defining boundary - * conditions of 1D grids. - * - * \param x Number of grid cells in x-direction - */ - BoundaryCondition(int x); - - /** - * Creates a new instance with 4 * max(x,y) elements. Used to describe the - * boundary conditions for 2D grids. NOTE: On non-squared grids there are more - * elements than needed and no exception is thrown if some index exceeding - * grid limits. - * - * QUESTION: why not use non-squared grids with the correct size? - * - * \param x Number of grid cells in x-direction - * \param y Number of grid cells in y-direction - * - */ - BoundaryCondition(int x, int y); - - /** - * Sets the boundary condition for a specific side of the grid. - * - * \param side Side for which the given boundary condition should be set. - * \param input_bc Instance of struct boundary_condition with desired boundary - * condition. - * - * \throws std::invalid_argument Indicates wrong dimensions of the grid. - * \throws std::out_of_range Indicates a out of range value for side. - */ - void setSide(uint8_t side, boundary_condition &input_bc); - - /** - * Sets the boundary condition for a specific side of the grid. - * - * \param side Side for which the given boundary condition should be set. - * \param input_bc Vector of boundary conditions for specific side. - * - * \throws std::invalid_argument Indicates wrong dimensions of the grid. - * \throws std::out_of_range Indicates a out of range value for side or - * invalid size of input vector. - */ - void setSide(uint8_t side, std::vector &input_bc); - - /** - * Returns a vector of boundary conditions of given side. Can be used to set - * custom boundary conditions and set back via setSide() with vector input. - * - * \param side Side which boundary conditions should be returned - * - * \returns Vector of boundary conditions - * - * \throws std::invalid_argument If given dimension is less or equal to 1. - * \throws std::out_of_range Indicates a out of range value for side. - */ - auto getSide(uint8_t side) -> std::vector; - - /** - * Get both boundary conditions of a given row (left and right). - * - * \param i Index of row - * - * \returns Left and right boundary values as an array (defined as data - * type bc_tuple). - */ - auto row_boundary(uint32_t i) const -> bc_tuple; - - /** - * Get both boundary conditions of a given column (top and bottom). - * - * \param i Index of column - * - * \returns Top and bottom boundary values as an array (defined as data - * type bc_tuple). - */ - auto col_boundary(uint32_t i) const -> bc_tuple; - - /** - * Sets one cell of the inner grid to a given boundary condition. - * - * \param bc Boundary condition to be set. - * \param x Index of grid cell in x direction. - * \param y Index of grid cell in y direction. - */ - void setInnerBC(boundary_condition bc, int x, int y); - - /** - * Unsets a previously set inner boundary condition. - * - * \param x Index of grid cell in x direction. - * \param y Index of grid cell in y direction. - */ - void unsetInnerBC(int x, int y); - - /** - * Returns the current boundary condition set for specific inner grid cell. - * - * \param x Index of grid cell in x direction. - * \param y Index of grid cell in y direction. - */ - auto getInnerBC(int x, int y) -> boundary_condition; - - /** - * Get a row of field and its inner boundary conditions. - * - * \param i Index of the row starting at 0. - * - * \returns Row of the inner boundary conditions of the field. - */ - auto getInnerRow(uint32_t i) const -> bc_vec; - - /** - * Get a column of field and its inner boundary conditions. - * - * \param i Index of the column starting at 0. - * - * \returns Column of the inner boundary conditions of the field. - */ - auto getInnerCol(uint32_t i) const -> bc_vec; - - /** - * Create an instance of boundary_condition data type. Can be seen as a helper - * function. - * - * \param type Type of the boundary condition. - * \param value According value of condition. - * - * \returns Instance of boundary_condition - */ - static boundary_condition returnBoundaryCondition(bctype type, double value) { - return {type, value}; - } - -private: - std::vector bc_internal; - - std::map inner_cells; - - uint8_t dim; - - std::array sizes; - uint32_t maxsize; - uint32_t maxindex; - - enum { X_DIM, Y_DIM }; - -// TODO combine the 'public' blocks -public: - /** - * Returns the left/right boundary condition for 1D grid. - * - * \param side Side of the boundary condition to get. - * - * \returns Boundary condition - */ - boundary_condition operator()(uint8_t side) const { - if (dim != 1) { - throw std::invalid_argument( - "Only 1D grid support 1 parameter in operator"); - } - if (side > 1) { - throw std::out_of_range("1D index out of range"); - } - return bc_internal[side]; - } - - /** - * Returns the boundary condition of a given side for 2D grids. - * - * \param side Side of the boundary condition to get. - * \param i Index of the boundary condition. - * - * \returns Boundary condition - */ - boundary_condition operator()(uint8_t side, uint32_t i) const { - if (dim != 2) { - throw std::invalid_argument( - "Only 2D grids support 2 parameters in operator"); - } - if (side > 3) { - throw std::out_of_range("2D index out of range"); - } - return bc_internal[side * maxsize + i]; - } - - /** - * Returns the left/right boundary condition for 1D grid. - * - * \param side Side of the boundary condition to get. - * - * \returns Boundary condition - */ - boundary_condition &operator()(uint8_t side) { - if (dim != 1) { - throw std::invalid_argument( - "Only 1D grid support 1 parameter in operator"); - } - if (side > 1) { - throw std::out_of_range("1D index out of range"); - } - return bc_internal[side]; - } - - /** - * Returns the boundary condition of a given side for 2D grids. - * - * \param side Side of the boundary condition to get. - * \param i Index of the boundary condition. - * - * \returns Boundary condition - */ - boundary_condition &operator()(uint8_t side, uint32_t i) { - if (dim != 2) { - throw std::invalid_argument("Explicit setting of bc value with 2 " - "parameters is only supported for 2D grids"); - } - if (side > 3) { - throw std::out_of_range("2D index out of range"); - } - return bc_internal[side * maxsize + i]; - } -}; - -} // namespace bc -} // namespace tug -#endif // BOUNDARYCONDITION_H_ diff --git a/include/tug/Diffusion.hpp b/include/tug/Diffusion.hpp deleted file mode 100644 index 52d4920..0000000 --- a/include/tug/Diffusion.hpp +++ /dev/null @@ -1,140 +0,0 @@ -#ifndef DIFFUSION_H_ -#define DIFFUSION_H_ - -#include "BoundaryCondition.hpp" -#include "Solver.hpp" -#include -#include -#include -#include - -namespace tug { -namespace diffusion { - -constexpr uint8_t MAX_ARR_SIZE = 3; - -/** - * Defines grid dimensions and boundary conditions. - * QUESTION: why is TugGrid a separate struct? For conceptual reasons? - * why not combine TugGrid and TugInput? - */ -typedef struct tug_grid_s { - std::array - grid_cells; /**< Count of grid cells in each of the 3 directions.*/ - std::array - domain_size; /**< Domain sizes in each of the 3 directions.*/ - bc::BoundaryCondition *bc = NULL; /**< Boundary conditions for the grid.*/ -} TugGrid; - -/** - * Besides containing the grid structure it holds also information about the - * desired time step to simulate and which solver to use. - */ -typedef struct tug_input_s { - double time_step; /**< Time step which should be simulated by diffusion.*/ - Eigen::VectorXd (*solver)(const Eigen::SparseMatrix &, - const Eigen::VectorXd &) = - tug::solver::ThomasAlgorithm; /**< Solver function to use.*/ - TugGrid grid; /**< Grid specification.*/ - - /** - * Set the desired time step for diffusion simulation. - * - * \param dt Time step in seconds. - */ - void setTimestep(double dt) { time_step = dt; } - - /** - * Set the count of grid cells in each dimension. - * - * \param x Count of grid cells in x direction. - * \param y Count of grid cells in y direction. - * \param z Count of grid cells in z direction. - */ - void setGridCellN(uint32_t x, uint32_t y = 0, uint32_t z = 0) { - grid.grid_cells[0] = x; - grid.grid_cells[1] = y; - grid.grid_cells[2] = z; - } - - /** - * Set the domain size of the grid in each direction. - - * \param Domain size in x direction. - * \param Domain size in y direction. - * \param Domain size in z direction. - */ - void setDomainSize(double x, double y = 0, double z = 0) { - grid.domain_size[0] = x; - grid.domain_size[1] = y; - grid.domain_size[2] = z; - } - - /** - * Set boundary conditions for grid instance. - * - * \param bc Boundary conditions to be set. - */ - void setBoundaryCondition(bc::BoundaryCondition &bc) { grid.bc = &bc; } - - /** - * Retrieve the set boundary condition from grid instance. - * - * \return Boundary condition object if boundary conditions were set, - * otherwise NULL. - */ - auto getBoundaryCondition() -> bc::BoundaryCondition { return *(grid.bc); } - - /** - * Set the solver function. - * - * \param f_in Pointer to function which takes a sparse matrix and a vector as - * input and returns another vector. - */ - void - setSolverFunction(Eigen::VectorXd (*f_in)(const Eigen::SparseMatrix &, - const Eigen::VectorXd &)) { - solver = f_in; - } -} TugInput; - -/** - * Solving 1D diffusion problems with Backward Time Centred Space scheme. - * - * \param input_param Object with information for the simulation e.g. grid - * dimensions, time step ... - * - * \param field Pointer to continious memory holding the concentrations for each - * grid cell of the grid. It doesn't matter if stored in row (most likely) or - * column major. - * - * \param alpha Pointer to continious memory holding the alpha for each grid - * cell of the grid. (NOTE: only constant alphas are supported currently) - * - * \return Runtime in seconds - */ -auto BTCS_1D(const TugInput &input_param, double *field, const double *alpha) - -> double; - -/** - * Solving 2D diffusion problems with Alternating-direction implicit method. - * - * \param input_param Object with information for the simulation e.g. grid - * dimensions, time step ... - * - * \param field Pointer to continious memory holding the concentrations for each - * grid cell of the grid. It doesn't matter if stored in row (most likely) or - * column major. - * - * \param alpha Pointer to continious memory holding the alpha for each grid - * cell of the grid. (NOTE: only constant alphas are supported currently) - * - * \return Runtime in seconds - */ -auto ADI_2D(const TugInput &input_param, double *field, const double *alpha) - -> double; - -} // namespace diffusion -} // namespace tug - -#endif // DIFFUSION_H_ diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index fb9f734..7bca259 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -16,7 +16,16 @@ using namespace std; */ enum APPROACH { FTCS_APPROACH, // Forward Time-Centered Space - BTCS_APPROACH // Backward Time-Centered Space + BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver +}; + +/** + * @brief Enum defining the Linear Equation solvers + * + */ +enum SOLVER { + EIGEN_LU_SOLVER, // EigenLU solver + THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for tridiagonal matrices }; /** @@ -124,6 +133,15 @@ class Simulation { */ void setIterations(int iterations); + /** + * @brief Set the desired linear equation solver to be used for BTCS approach. Without effect + * in case of FTCS approach. + * + * @param solver Solver to be used. Default is Thomas Algorithm as it is more efficient for + * tridiagonal Matrices. + */ + void setSolver(SOLVER solver); + /** * @brief Return the currently set iterations to be calculated. * @@ -174,5 +192,6 @@ class Simulation { Grid &grid; Boundary &bc; APPROACH approach; + SOLVER solver; }; diff --git a/include/tug/Solver.hpp b/include/tug/Solver.hpp deleted file mode 100644 index a7abd9d..0000000 --- a/include/tug/Solver.hpp +++ /dev/null @@ -1,41 +0,0 @@ -#ifndef SOLVER_H_ -#define SOLVER_H_ - -#include -#include - -namespace tug { -namespace solver { - -/** - * Solving linear equation system with LU-decomposition implemented in Eigen - * library. - * - * \param A_matrix The A matrix represented as a sparse matrix using Eigen - * library. - * - * \param b_vector Right hand side vector of the linear equation system. - * - * \return Solution represented as vector. - */ -auto EigenLU(const Eigen::SparseMatrix &A_matrix, - const Eigen::VectorXd &b_vector) -> Eigen::VectorXd; - -/** - * Solving linear equation system with brutal implementation of the Thomas - * algorithm (a.k.a. Tridiagonal matrix algorithm). - * - * \param A_matrix The A matrix represented as a sparse matrix using Eigen - * library. - * - * \param b_vector Right hand side vector of the linear equation system. - * - * \return Solution represented as vector. - */ -auto ThomasAlgorithm(const Eigen::SparseMatrix &A_matrix, - const Eigen::VectorXd &b_vector) -> Eigen::VectorXd; - -} // namespace solver -} // namespace tug - -#endif // SOLVER_H_ diff --git a/include/tug/progressbar.hpp b/include/tug/progressbar.hpp deleted file mode 100644 index 50ab6d0..0000000 --- a/include/tug/progressbar.hpp +++ /dev/null @@ -1,191 +0,0 @@ -// The MIT License (MIT) -// -// Copyright (c) 2019 Luigi Pertoldi -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to -// deal in the Software without restriction, including without limitation the -// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or -// sell copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in -// all copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS -// IN THE SOFTWARE. -// -// ============================================================================ -// ___ ___ ___ __ ___ ____ __ __ ___ __ ___ -// | |_) | |_) / / \ / /`_ | |_) | |_ ( (` ( (` | |_) / /\ | |_) -// |_| |_| \ \_\_/ \_\_/ |_| \ |_|__ _)_) _)_) |_|_) /_/--\ |_| \_ -// -// Very simple progress bar for c++ loops with internal running variable -// -// Author: Luigi Pertoldi -// Created: 3 dic 2016 -// -// Notes: The bar must be used when there's no other possible source of output -// inside the for loop -// - -#ifndef __PROGRESSBAR_HPP -#define __PROGRESSBAR_HPP - -#include -#include -#include -#include - -class progressbar { - - public: - // default destructor - ~progressbar() = default; - - // delete everything else - progressbar (progressbar const&) = delete; - progressbar& operator=(progressbar const&) = delete; - progressbar (progressbar&&) = delete; - progressbar& operator=(progressbar&&) = delete; - - // default constructor, must call set_niter later - inline progressbar(); - inline progressbar(int n, bool showbar=true, std::ostream& out=std::cerr); - - // reset bar to use it again - inline void reset(); - // set number of loop iterations - inline void set_niter(int iter); - // chose your style - inline void set_done_char(const std::string& sym) {done_char = sym;} - inline void set_todo_char(const std::string& sym) {todo_char = sym;} - inline void set_opening_bracket_char(const std::string& sym) {opening_bracket_char = sym;} - inline void set_closing_bracket_char(const std::string& sym) {closing_bracket_char = sym;} - // to show only the percentage - inline void show_bar(bool flag = true) {do_show_bar = flag;} - // set the output stream - inline void set_output_stream(const std::ostream& stream) {output.rdbuf(stream.rdbuf());} - // main function - inline void update(); - - private: - int progress; - int n_cycles; - int last_perc; - bool do_show_bar; - bool update_is_called; - - std::string done_char; - std::string todo_char; - std::string opening_bracket_char; - std::string closing_bracket_char; - - std::ostream& output; -}; - -inline progressbar::progressbar() : - progress(0), - n_cycles(0), - last_perc(0), - do_show_bar(true), - update_is_called(false), - done_char("#"), - todo_char(" "), - opening_bracket_char("["), - closing_bracket_char("]"), - output(std::cerr) {} - -inline progressbar::progressbar(int n, bool showbar, std::ostream& out) : - progress(0), - n_cycles(n), - last_perc(0), - do_show_bar(showbar), - update_is_called(false), - done_char("#"), - todo_char(" "), - opening_bracket_char("["), - closing_bracket_char("]"), - output(out) {} - -inline void progressbar::reset() { - progress = 0, - update_is_called = false; - last_perc = 0; - return; -} - -inline void progressbar::set_niter(int niter) { - if (niter <= 0) throw std::invalid_argument( - "progressbar::set_niter: number of iterations null or negative"); - n_cycles = niter; - return; -} - -inline void progressbar::update() { - - if (n_cycles == 0) throw std::runtime_error( - "progressbar::update: number of cycles not set"); - - if (!update_is_called) { - if (do_show_bar == true) { - output << opening_bracket_char; - for (int _ = 0; _ < 50; _++) output << todo_char; - output << closing_bracket_char << " 0%"; - } - else output << "0%"; - } - update_is_called = true; - - int perc = 0; - - // compute percentage, if did not change, do nothing and return - perc = progress*100./(n_cycles-1); - if (perc < last_perc) return; - - // update percentage each unit - if (perc == last_perc + 1) { - // erase the correct number of characters - if (perc <= 10) output << "\b\b" << perc << '%'; - else if (perc > 10 and perc < 100) output << "\b\b\b" << perc << '%'; - else if (perc == 100) output << "\b\b\b" << perc << '%'; - } - if (do_show_bar == true) { - // update bar every ten units - if (perc % 2 == 0) { - // erase closing bracket - output << std::string(closing_bracket_char.size(), '\b'); - // erase trailing percentage characters - if (perc < 10) output << "\b\b\b"; - else if (perc >= 10 && perc < 100) output << "\b\b\b\b"; - else if (perc == 100) output << "\b\b\b\b\b"; - - // erase 'todo_char' - for (int j = 0; j < 50-(perc-1)/2; ++j) { - output << std::string(todo_char.size(), '\b'); - } - - // add one additional 'done_char' - if (perc == 0) output << todo_char; - else output << done_char; - - // refill with 'todo_char' - for (int j = 0; j < 50-(perc-1)/2-1; ++j) output << todo_char; - - // readd trailing percentage characters - output << closing_bracket_char << ' ' << perc << '%'; - } - } - last_perc = perc; - ++progress; - output << std::flush; - - return; -} - -#endif diff --git a/src/BTCS.cpp b/src/BTCS.cpp deleted file mode 100644 index 9e31d52..0000000 --- a/src/BTCS.cpp +++ /dev/null @@ -1,328 +0,0 @@ -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#include "TugUtils.hpp" - -#ifdef _OPENMP -#include -#else -#define omp_get_thread_num() 0 -#endif - -inline auto -init_delta(const std::array &domain_size, - const std::array &grid_cells, - const uint8_t dim) -> std::vector { - std::vector out(dim); - for (uint8_t i = 0; i < dim; i++) { - // calculate 'size' of each cell in grid - out[i] = (double)(domain_size.at(i) / grid_cells.at(i)); - } - return out; -} - -namespace { -enum { GRID_1D = 1, GRID_2D, GRID_3D }; - -constexpr int BTCS_MAX_DEP_PER_CELL = 3; -constexpr int BTCS_2D_DT_SIZE = 2; - -using DMatrixRowMajor = - Eigen::Matrix; -using DVectorRowMajor = - Eigen::Matrix; - -inline auto getBCFromFlux(tug::bc::boundary_condition bc, double neighbor_c, - double neighbor_alpha) -> double { - double val = 0; - - if (bc.type == tug::bc::BC_TYPE_CLOSED) { - val = neighbor_c; - } else if (bc.type == tug::bc::BC_TYPE_FLUX) { - // TODO - // val = bc[index].value; - } else { - // TODO: implement error handling here. Type was set to wrong value. - } - - return val; -} - -auto calc_d_ortho(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha, - const tug::bc::BoundaryCondition &bc, bool transposed, - double time_step, double dx) -> DMatrixRowMajor { - - uint8_t upper = (transposed ? tug::bc::BC_SIDE_LEFT : tug::bc::BC_SIDE_TOP); - uint8_t lower = - (transposed ? tug::bc::BC_SIDE_RIGHT : tug::bc::BC_SIDE_BOTTOM); - - int n_rows = c.rows(); - int n_cols = c.cols(); - - DMatrixRowMajor d_ortho(n_rows, n_cols); - - std::array y_values{}; - - // first, iterate over first row - for (int j = 0; j < n_cols; j++) { - tug::bc::boundary_condition tmp_bc = bc(upper, j); - double sy = (time_step * alpha(0, j)) / (dx * dx); - - y_values[0] = (tmp_bc.type == tug::bc::BC_TYPE_CONSTANT - ? tmp_bc.value - : getBCFromFlux(tmp_bc, c(0, j), alpha(0, j))); - y_values[1] = c(0, j); - y_values[2] = c(1, j); - - d_ortho(0, j) = -sy * (2 * y_values[0] - 3 * y_values[1] + y_values[2]); - } - - // then iterate over inlet - for (int i = 1; i < n_rows - 1; i++) { - for (int j = 0; j < n_cols; j++) { - double sy = (time_step * alpha(i, j)) / (dx * dx); - - y_values[0] = c(i - 1, j); - y_values[1] = c(i, j); - y_values[2] = c(i + 1, j); - - d_ortho(i, j) = -sy * (y_values[0] - 2 * y_values[1] + y_values[2]); - } - } - - int end = n_rows - 1; - - // and finally over last row - for (int j = 0; j < n_cols; j++) { - tug::bc::boundary_condition tmp_bc = bc(lower, j); - double sy = (time_step * alpha(end, j)) / (dx * dx); - - y_values[0] = c(end - 1, j); - y_values[1] = c(end, j); - y_values[2] = (tmp_bc.type == tug::bc::BC_TYPE_CONSTANT - ? tmp_bc.value - : getBCFromFlux(tmp_bc, c(end, j), alpha(end, j))); - - d_ortho(end, j) = -sy * (y_values[0] - 3 * y_values[1] + 2 * y_values[2]); - } - - return d_ortho; -} - -auto fillMatrixFromRow(const DVectorRowMajor &alpha, - const tug::bc::bc_vec &bc_inner, int size, double dx, - double time_step) -> Eigen::SparseMatrix { - - Eigen::SparseMatrix A_matrix(size + 2, size + 2); - - A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, BTCS_MAX_DEP_PER_CELL)); - - double sx = 0; - - int A_size = A_matrix.cols(); - - A_matrix.insert(0, 0) = 1; - - if (bc_inner[0].type != tug::bc::BC_UNSET) { - if (bc_inner[0].type != tug::bc::BC_TYPE_CONSTANT) { - throw_invalid_argument("Inner boundary conditions with other type than " - "BC_TYPE_CONSTANT are currently not supported."); - } - A_matrix.insert(1, 1) = 1; - } else { // if BC_UNSET - sx = (alpha[0] * time_step) / (dx * dx); - A_matrix.insert(1, 1) = -1. - 3. * sx; - A_matrix.insert(1, 0) = 2. * sx; - A_matrix.insert(1, 2) = sx; - } - - for (int j = 2, k = j - 1; k < size - 1; j++, k++) { - if (bc_inner[k].type != tug::bc::BC_UNSET) { - if (bc_inner[k].type != tug::bc::BC_TYPE_CONSTANT) { - throw_invalid_argument("Inner boundary conditions with other type than " - "BC_TYPE_CONSTANT are currently not supported."); - } - A_matrix.insert(j, j) = 1; - continue; - } // if BC_UNSET - sx = (alpha[k] * time_step) / (dx * dx); - - A_matrix.insert(j, j) = -1. - 2. * sx; - A_matrix.insert(j, (j - 1)) = sx; - A_matrix.insert(j, (j + 1)) = sx; - } - - if (bc_inner[size - 1].type != tug::bc::BC_UNSET) { - if (bc_inner[size - 1].type != tug::bc::BC_TYPE_CONSTANT) { - throw_invalid_argument("Inner boundary conditions with other type than " - "BC_TYPE_CONSTANT are currently not supported."); - } - A_matrix.insert(A_size - 2, A_size - 2) = 1; - } else { // if BC_UNSET - sx = (alpha[size - 1] * time_step) / (dx * dx); - A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx; - A_matrix.insert(A_size - 2, A_size - 3) = sx; - A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx; - } - - A_matrix.insert(A_size - 1, A_size - 1) = 1; - - return A_matrix; -} - -auto fillVectorFromRow(const DVectorRowMajor &c, const DVectorRowMajor &alpha, - const tug::bc::bc_tuple &bc, - const tug::bc::bc_vec &bc_inner, - const DVectorRowMajor &d_ortho, int size, double dx, - double time_step) -> Eigen::VectorXd { - - Eigen::VectorXd b_vector(size + 2); - - tug::bc::boundary_condition left = bc[0]; - tug::bc::boundary_condition right = bc[1]; - - bool left_constant = (left.type == tug::bc::BC_TYPE_CONSTANT); - bool right_constant = (right.type == tug::bc::BC_TYPE_CONSTANT); - - int b_size = b_vector.size(); - - for (int j = 0; j < size; j++) { - if (bc_inner[j].type != tug::bc::BC_UNSET) { - if (bc_inner[j].type != tug::bc::BC_TYPE_CONSTANT) { - throw_invalid_argument("Inner boundary conditions with other type than " - "BC_TYPE_CONSTANT are currently not supported."); - } - b_vector[j + 1] = bc_inner[j].value; - continue; - } // if BC_UNSET - b_vector[j + 1] = -c[j] + d_ortho[j]; - } - - // this is not correct currently.We will fix this when we are able to define - // FLUX boundary conditions - b_vector[0] = - (left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0])); - - b_vector[b_size - 1] = - (right_constant ? right.value - : getBCFromFlux(right, c[size - 1], alpha[size - 1])); - - return b_vector; -} - -auto setupBTCSAndSolve( - DVectorRowMajor &c, const tug::bc::bc_tuple bc_ghost, - const tug::bc::bc_vec &bc_inner, const DVectorRowMajor &alpha, double dx, - double time_step, int size, const DVectorRowMajor &d_ortho, - Eigen::VectorXd (*solver)(const Eigen::SparseMatrix &, - const Eigen::VectorXd &)) -> DVectorRowMajor { - - const Eigen::SparseMatrix A_matrix = - fillMatrixFromRow(alpha, bc_inner, size, dx, time_step); - - const Eigen::VectorXd b_vector = fillVectorFromRow( - c, alpha, bc_ghost, bc_inner, d_ortho, size, dx, time_step); - - // solving of the LEQ - Eigen::VectorXd x_vector = solver(A_matrix, b_vector); - - DVectorRowMajor out_vector = x_vector.segment(1, size); - - return out_vector; -} - -} // namespace - // -auto tug::diffusion::BTCS_1D(const tug::diffusion::TugInput &input_param, - double *field, const double *alpha) -> double { - - auto start = time_marker(); - - uint32_t size = input_param.grid.grid_cells[0]; - - auto deltas = init_delta(input_param.grid.domain_size, - input_param.grid.grid_cells, GRID_1D); - double dx = deltas[0]; - - double time_step = input_param.time_step; - - const tug::bc::BoundaryCondition bc = - (input_param.grid.bc != nullptr ? *input_param.grid.bc - : tug::bc::BoundaryCondition(size)); - - Eigen::Map c_in(field, size); - Eigen::Map alpha_in(alpha, size); - - DVectorRowMajor input_field = c_in.row(0); - - DVectorRowMajor output = setupBTCSAndSolve( - input_field, bc.row_boundary(0), bc.getInnerRow(0), alpha_in, dx, - time_step, size, Eigen::VectorXd::Constant(size, 0), input_param.solver); - - c_in.row(0) << output; - - auto end = time_marker(); - - return diff_time(start, end); -} - -auto tug::diffusion::ADI_2D(const tug::diffusion::TugInput &input_param, - double *field, const double *alpha) -> double { - - auto start = time_marker(); - - uint32_t n_cols = input_param.grid.grid_cells[0]; - uint32_t n_rows = input_param.grid.grid_cells[1]; - - auto deltas = init_delta(input_param.grid.domain_size, - input_param.grid.grid_cells, GRID_2D); - double dx = deltas[0]; - double dy = deltas[1]; - - double local_dt = input_param.time_step / BTCS_2D_DT_SIZE; - - tug::bc::BoundaryCondition bc = - (input_param.grid.bc != nullptr - ? *input_param.grid.bc - : tug::bc::BoundaryCondition(n_cols, n_rows)); - - Eigen::Map c_in(field, n_rows, n_cols); - Eigen::Map alpha_in(alpha, n_rows, n_cols); - - DMatrixRowMajor d_ortho = - calc_d_ortho(c_in, alpha_in, bc, false, local_dt, dx); - -#pragma omp parallel for schedule(dynamic) - for (int i = 0; i < n_rows; i++) { - DVectorRowMajor input_field = c_in.row(i); - DVectorRowMajor output = setupBTCSAndSolve( - input_field, bc.row_boundary(i), bc.getInnerRow(i), alpha_in.row(i), dx, - local_dt, n_cols, d_ortho.row(i), input_param.solver); - c_in.row(i) << output; - } - - d_ortho = calc_d_ortho(c_in.transpose(), alpha_in.transpose(), bc, true, - local_dt, dy); - -#pragma omp parallel for schedule(dynamic) - for (int i = 0; i < n_cols; i++) { - DVectorRowMajor input_field = c_in.col(i); - DVectorRowMajor output = setupBTCSAndSolve( - input_field, bc.col_boundary(i), bc.getInnerCol(i), alpha_in.col(i), dy, - local_dt, n_rows, d_ortho.row(i), input_param.solver); - c_in.col(i) << output.transpose(); - } - - auto end = time_marker(); - - return diff_time(start, end); -} diff --git a/src/BTCSv2.cpp b/src/BTCSv2.cpp index f8e6ea6..8751a43 100644 --- a/src/BTCSv2.cpp +++ b/src/BTCSv2.cpp @@ -264,7 +264,8 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, // solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector -static VectorXd solve(SparseMatrix &A, VectorXd &b) { +// use of EigenLU solver +static VectorXd EigenLUAlgorithm(SparseMatrix &A, VectorXd &b) { SparseLU> solver; solver.analyzePattern(A); @@ -273,9 +274,53 @@ static VectorXd solve(SparseMatrix &A, VectorXd &b) { return solver.solve(b); } +// solver for linear equation system; A corresponds to coefficient matrix, +// b to the solution vector +// implementation of Thomas Algorithm +static VectorXd ThomasAlgorithm(SparseMatrix &A, VectorXd &b) { + uint32_t n = b.size(); + + Eigen::VectorXd a_diag(n); + Eigen::VectorXd b_diag(n); + Eigen::VectorXd c_diag(n); + Eigen::VectorXd x_vec = b; + + // Fill diagonals vectors + b_diag[0] = A.coeff(0, 0); + c_diag[0] = A.coeff(0, 1); + + for (int i = 1; i < n - 1; i++) { + a_diag[i] = A.coeff(i, i - 1); + b_diag[i] = A.coeff(i, i); + c_diag[i] = A.coeff(i, i + 1); + } + a_diag[n - 1] = A.coeff(n - 1, n - 2); + b_diag[n - 1] = A.coeff(n - 1, n - 1); + + // start solving - c_diag and x_vec are overwritten + n--; + c_diag[0] /= b_diag[0]; + x_vec[0] /= b_diag[0]; + + for (int i = 1; i < n; i++) { + c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i - 1]; + x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i - 1]) / + (b_diag[i] - a_diag[i] * c_diag[i - 1]); + } + + x_vec[n] = (x_vec[n] - a_diag[n] * x_vec[n - 1]) / + (b_diag[n] - a_diag[n] * c_diag[n - 1]); + + for (int i = n; i-- > 0;) { + x_vec[i] -= c_diag[i] * x_vec[i + 1]; + } + + return x_vec; +} + // BTCS solution for 1D grid -static void BTCS_1D(Grid &grid, Boundary &bc, double ×tep) { +static void BTCS_1D(Grid &grid, Boundary &bc, double ×tep, VectorXd (*solverFunc) (SparseMatrix &A, VectorXd &b)) { int length = grid.getLength(); double sx = timestep / (grid.getDelta() * grid.getDelta()); @@ -301,7 +346,7 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double ×tep) { b(length-1) += 2 * sx * alpha(0,length-1) * bcRight[0].getValue(); } - concentrations_t1 = solve(A, b); + concentrations_t1 = solverFunc(A, b); for (int j = 0; j < length; j++) { concentrations(0,j) = concentrations_t1(j); @@ -312,7 +357,7 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double ×tep) { // BTCS solution for 2D grid -static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep) { +static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep, VectorXd (*solverFunc) (SparseMatrix &A, VectorXd &b)) { int rowMax = grid.getRow(); int colMax = grid.getCol(); double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); @@ -342,7 +387,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep) { SparseLU> solver; - row_t1 = solve(A, b); + row_t1 = solverFunc(A, b); concentrations_t1.row(i) = row_t1; } @@ -360,7 +405,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep) { b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, bcLeft, bcRight, rowMax, i, sy, sx); - row_t1 = solve(A, b); + row_t1 = solverFunc(A, b); concentrations.row(i) = row_t1; } @@ -371,12 +416,23 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep) { } -// entry point; differentiate between 1D and 2D grid -static void BTCS(Grid &grid, Boundary &bc, double ×tep) { +// entry point for EigenLU solver; differentiate between 1D and 2D grid +static void BTCS_LU(Grid &grid, Boundary &bc, double ×tep) { if (grid.getDim() == 1) { - BTCS_1D(grid, bc, timestep); + BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); } else if (grid.getDim() == 2) { - BTCS_2D(grid, bc, timestep); + BTCS_2D(grid, bc, timestep, EigenLUAlgorithm); + } else { + throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); + } +} + +// entry point for Thomas algorithm solver; differentiate 1D and 2D grid +static void BTCS_Thomas(Grid &grid, Boundary &bc, double ×tep) { + if (grid.getDim() == 1) { + BTCS_1D(grid, bc, timestep, ThomasAlgorithm); + } else if (grid.getDim() == 2) { + BTCS_2D(grid, bc, timestep, ThomasAlgorithm); } else { throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); } diff --git a/src/Boundary.cpp b/src/Boundary.cpp index bba7d9d..4f9890b 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -1,5 +1,4 @@ #include "TugUtils.hpp" -#include "tug/BoundaryCondition.hpp" #include #include #include @@ -129,7 +128,7 @@ VectorXd Boundary::getBoundarySideValues(BC_SIDE side) { VectorXd values(length); for (int i = 0; i < length; i++) { - if (getBoundaryElementType(side, i) == tug::bc::BC_TYPE_CLOSED) { + if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { values(i) = -1; continue; } diff --git a/src/BoundaryCondition.cpp b/src/BoundaryCondition.cpp deleted file mode 100644 index 63f5bc8..0000000 --- a/src/BoundaryCondition.cpp +++ /dev/null @@ -1,209 +0,0 @@ -#include -#include -#include - -#include "TugUtils.hpp" - -constexpr uint8_t DIM_1D = 2; -constexpr uint8_t DIM_2D = 4; - -tug::bc::BoundaryCondition::BoundaryCondition(int x) { - this->bc_internal.resize(DIM_1D, {0, 0}); - this->dim = 1; - // this value is actually unused - this->maxsize = 1; - - this->sizes[X_DIM] = x; - this->sizes[Y_DIM] = 1; - - this->maxindex = x - 1; -} - -tug::bc::BoundaryCondition::BoundaryCondition(int x, int y) { - this->maxsize = (x >= y ? x : y); - this->bc_internal.resize(DIM_2D * maxsize, {0, 0}); - this->dim = 2; - - this->sizes[X_DIM] = x; - this->sizes[Y_DIM] = y; - - this->maxindex = (x * y) - 1; -} - -void tug::bc::BoundaryCondition::setSide( - uint8_t side, tug::bc::boundary_condition &input_bc) { - // QUESTION: why cant the BC be changed for dim = 1? - if (this->dim == 1) { - throw_invalid_argument("setSide requires at least a 2D grid"); - } - if (side > 3) { - throw_out_of_range("Invalid range for 2D grid"); - } - - uint32_t size = - (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT - ? this->sizes[Y_DIM] - : this->sizes[X_DIM]); - - for (uint32_t i = 0; i < size; i++) { - this->bc_internal[side * maxsize + i] = input_bc; - } -} - -void tug::bc::BoundaryCondition::setSide( - uint8_t side, std::vector &input_bc) { - if (this->dim == 1) { - throw_invalid_argument("setSide requires at least a 2D grid"); - } - if (side > 3) { - throw_out_of_range("Invalid range for 2D grid"); - } - - uint32_t size = - (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT - ? this->sizes[Y_DIM] - : this->sizes[X_DIM]); - - if (input_bc.size() > size) { - throw_out_of_range("Input vector is greater than maximum excpected value"); - } - - for (int i = 0; i < size; i++) { - bc_internal[this->maxsize * side + i] = input_bc[i]; - } -} - -auto tug::bc::BoundaryCondition::getSide(uint8_t side) - -> std::vector { - if (this->dim == 1) { - throw_invalid_argument("getSide requires at least a 2D grid"); - } - if (side > 3) { - throw_out_of_range("Invalid range for 2D grid"); - } - - uint32_t size = - (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT - ? this->sizes[Y_DIM] - : this->sizes[X_DIM]); - - std::vector out(size); - - for (int i = 0; i < size; i++) { - out[i] = this->bc_internal[this->maxsize * side + i]; - } - - return out; -} - -auto tug::bc::BoundaryCondition::row_boundary(uint32_t i) const - -> tug::bc::bc_tuple { - if (i >= this->sizes[Y_DIM]) { - throw_out_of_range("Index out of range"); - } - - return {this->bc_internal[BC_SIDE_LEFT * this->maxsize + i], - this->bc_internal[BC_SIDE_RIGHT * this->maxsize + i]}; -} - -auto tug::bc::BoundaryCondition::col_boundary(uint32_t i) const - -> tug::bc::bc_tuple { - if (this->dim == 1) { - throw_invalid_argument("Access of column requires at least 2D grid"); - } - if (i >= this->sizes[X_DIM]) { - throw_out_of_range("Index out of range"); - } - - return {this->bc_internal[BC_SIDE_TOP * this->maxsize + i], - this->bc_internal[BC_SIDE_BOTTOM * this->maxsize + i]}; -} - -auto tug::bc::BoundaryCondition::getInnerRow(uint32_t i) const -> bc_vec { - if (i >= this->sizes[Y_DIM]) { - throw_out_of_range("Index is out of range"); - } - - bc_vec row(this->sizes[X_DIM], {tug::bc::BC_UNSET, 0}); - - if (this->inner_cells.empty()) { - return row; - } - - uint32_t index_min = i * this->sizes[X_DIM]; - uint32_t index_max = ((i + 1) * this->sizes[X_DIM]) - 1; - - for (auto const &cell : this->inner_cells) { - if (cell.first < index_min) { - continue; - } - if (cell.first > index_max) { - break; - } - - row[cell.first - index_min] = cell.second; - } - - return row; -} - -auto tug::bc::BoundaryCondition::getInnerCol(uint32_t i) const -> bc_vec { - if (this->dim != 2) { - throw_invalid_argument("getInnerCol is only applicable for 2D grids"); - } - if (i >= this->sizes[X_DIM]) { - throw_out_of_range("Index is out of range"); - } - - bc_vec col(this->sizes[Y_DIM], {tug::bc::BC_UNSET, 0}); - - if (this->inner_cells.empty()) { - return col; - } - - for (auto const &cell : this->inner_cells) { - if (cell.first % this->sizes[X_DIM] == i) { - col[cell.first / this->sizes[X_DIM]] = cell.second; - } - } - - return col; -} - -void tug::bc::BoundaryCondition::setInnerBC(boundary_condition bc, int x, - int y = 0) { - if (x >= this->sizes[X_DIM] || y >= this->sizes[Y_DIM]) { - throw_out_of_range("One input parameter is out of range"); - } - uint32_t index = y * this->sizes[X_DIM] + x; - auto it = this->inner_cells.find(index); - - if (it != this->inner_cells.end()) { - it->second = bc; - return; - } - - this->inner_cells.insert({index, bc}); -} - -void tug::bc::BoundaryCondition::unsetInnerBC(int x, int y) { - uint32_t index = y * this->sizes[X_DIM] + x; - this->inner_cells.erase(index); -} - -auto tug::bc::BoundaryCondition::getInnerBC(int x, int y = 0) - -> boundary_condition { - if (x >= this->sizes[X_DIM] || y >= this->sizes[Y_DIM]) { - throw_out_of_range("One input parameter is out of range"); - } - - uint32_t index = y * this->sizes[X_DIM] + x; - - auto it = this->inner_cells.find(index); - - if (it != this->inner_cells.end()) { - return it->second; - } - - return {tug::bc::BC_UNSET, 0}; -} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5f1fca0..b7aceb2 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,4 @@ -add_library(tug BTCS.cpp BoundaryCondition.cpp Solver.cpp Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCSv2.cpp) +add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCSv2.cpp) target_link_libraries(tug Eigen3::Eigen) diff --git a/src/README.md b/src/README.md index d626308..82cd165 100644 --- a/src/README.md +++ b/src/README.md @@ -5,15 +5,12 @@ This is the src-directory that holds the source code to the implementation of th src/ ├── CMakeFiles/ ├── Boundary.cpp - ├── BoundaryCondition.cpp - ├── BTCS.cpp ├── BTCSv2.cpp ├── CMakeLists.txt ├── FTCS.cpp ├── Grid.cpp ├── README.md ├── Simulation.cpp - ├── Solver.cpp └── TugUtils.hpp @@ -23,10 +20,6 @@ src/ **Boundary.cpp** Implementation of Boundary class, that holds the boundary conditions. -**BoundaryCondition.cpp** Outdated implementation of boundary conditions. - -**BTCS.cpp** Outdated implementation of BTCS solution to homogeneous diffusion. - **BTCSv2.cpp** Implementation of BTCS solution to heterogeneous diffusion in 1D and 2D. **CMakeLists.txt** CMakeLists for this directory. @@ -39,7 +32,5 @@ src/ **Simulation.cpp** Implementation of Simulation class, that holds all of the information for a specific simulation run, as well as the Boundary and Grid objects. -**Solver.cpp** Outdated implementation of Eigen solvers. - **TugUtils.hpp** Helper functions for other source files. diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 17645a5..cbdfee9 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -7,8 +7,10 @@ #include +#ifndef SIMULATION_H_ +#define SIMULATION_H_ + #include "BTCSv2.cpp" -#include using namespace std; @@ -16,6 +18,7 @@ using namespace std; Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid), bc(bc) { this->approach = approach; + this->solver = THOMAS_ALGORITHM_SOLVER; this->timestep = -1; // error per default this->iterations = -1; this->innerIterations = 1; @@ -54,53 +57,62 @@ void Simulation::setTimestep(double timestep) { throw_invalid_argument("Timestep has to be greater than zero."); } - double deltaRowSquare; - double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - double minDeltaSquare; - double maxAlphaX, maxAlphaY, maxAlpha; - if (grid.getDim() == 2) { + if (approach == FTCS_APPROACH) { - deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + double deltaRowSquare; + double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + double minDeltaSquare; + double maxAlphaX, maxAlphaY, maxAlpha; + string dim; + if (grid.getDim() == 2) { + dim = "2D"; - minDeltaSquare = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - maxAlphaX = grid.getAlphaX().maxCoeff(); - maxAlphaY = grid.getAlphaY().maxCoeff(); - maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; + deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); - } else if (grid.getDim() == 1) { - minDeltaSquare = deltaColSquare; - maxAlpha = grid.getAlpha().maxCoeff(); + minDeltaSquare = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; + maxAlphaX = grid.getAlphaX().maxCoeff(); + maxAlphaY = grid.getAlphaY().maxCoeff(); + maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; - - } else { - throw_invalid_argument("Critical error: Undefined number of dimensions!"); - } + } else if (grid.getDim() == 1) { + dim = "1D"; + minDeltaSquare = deltaColSquare; + maxAlpha = grid.getAlpha().maxCoeff(); + + } else { + throw_invalid_argument("Critical error: Undefined number of dimensions!"); + } - // TODO check formula 1D case - double CFL_MDL = minDeltaSquare / (4*maxAlpha); // Formula from Marco --> seems to be unstable - double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia + // Courant-Friedrichs-Lewy condition + double cfl = minDeltaSquare / (4*maxAlpha); - cout << "FTCS_2D :: CFL condition MDL: " << CFL_MDL << endl; - // cout << "FTCS_2D :: CFL condition Wiki: " << CFL_Wiki << endl; - cout << "FTCS_2D :: required dt=" << timestep << endl; + // stability equation from Wikipedia; might be useful if applied cfl does not work in some cases + // double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); - if (timestep > CFL_MDL) { + cout << "FTCS_" << dim << " :: CFL condition MDL: " << cfl << endl; + cout << "FTCS_" << dim << " :: required dt=" << timestep << endl; - this->innerIterations = (int)ceil(timestep / CFL_MDL); - this->timestep = timestep / (double)innerIterations; + if (timestep > cfl) { - cerr << "Warning: Timestep was adjusted, because of stability " - "conditions. Time duration was approximately preserved by " - "adjusting internal number of iterations." - << endl; - cout << "FTCS_2D :: Required " << this->innerIterations - << " inner iterations with dt=" << this->timestep << endl; + this->innerIterations = (int)ceil(timestep / cfl); + this->timestep = timestep / (double)innerIterations; + + cerr << "Warning :: Timestep was adjusted, because of stability " + "conditions. Time duration was approximately preserved by " + "adjusting internal number of iterations." + << endl; + cout << "FTCS_" << dim << " :: Required " << this->innerIterations + << " inner iterations with dt=" << this->timestep << endl; + + } else { + + this->timestep = timestep; + cout << "FTCS_" << dim << " :: No inner iterations required, dt=" << timestep << endl; + + } } else { - this->timestep = timestep; - cout << "FTCS_2D :: No inner iterations required, dt=" << timestep << endl; - } } @@ -116,6 +128,16 @@ void Simulation::setIterations(int iterations) { this->iterations = iterations; } +void Simulation::setSolver(SOLVER solver) { + if (this->approach == FTCS_APPROACH) { + cerr << "Warning: Solver was set, but FTCS approach initialized. Setting the solver " + "is thus without effect." + << endl; + } + + this->solver = solver; +} + int Simulation::getIterations() { return this->iterations; } @@ -196,7 +218,7 @@ void Simulation::run() { auto begin = std::chrono::high_resolution_clock::now(); - if (approach == FTCS_APPROACH) { + if (approach == FTCS_APPROACH) { // FTCS case for (int i = 0; i < iterations * innerIterations; i++) { if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { printConcentrationsConsole(); @@ -215,18 +237,32 @@ void Simulation::run() { } } - } else if (approach == BTCS_APPROACH) { + } else if (approach == BTCS_APPROACH) { // BTCS case - for (int i = 0; i < iterations * innerIterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } + if (solver == EIGEN_LU_SOLVER) { + for (int i = 0; i < iterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } - BTCS(this->grid, this->bc, this->timestep); - + BTCS_LU(this->grid, this->bc, this->timestep); + + } + } else if (solver == THOMAS_ALGORITHM_SOLVER) { + for (int i = 0; i < iterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + BTCS_Thomas(this->grid, this->bc, this->timestep); + + } } } @@ -247,3 +283,5 @@ void Simulation::run() { } } + +#endif diff --git a/src/Solver.cpp b/src/Solver.cpp deleted file mode 100644 index 4733c8f..0000000 --- a/src/Solver.cpp +++ /dev/null @@ -1,61 +0,0 @@ -#include - -#include -#include - -auto tug::solver::EigenLU(const Eigen::SparseMatrix &A_matrix, - const Eigen::VectorXd &b_vector) -> Eigen::VectorXd { - - Eigen::SparseLU, Eigen::COLAMDOrdering> - solver; - solver.analyzePattern(A_matrix); - - solver.factorize(A_matrix); - - Eigen::VectorXd x_vec = solver.solve(b_vector); - - return solver.solve(b_vector); -} - -auto tug::solver::ThomasAlgorithm(const Eigen::SparseMatrix &A_matrix, - const Eigen::VectorXd &b_vector) - -> Eigen::VectorXd { - uint32_t n = b_vector.size(); - - Eigen::VectorXd a_diag(n); - Eigen::VectorXd b_diag(n); - Eigen::VectorXd c_diag(n); - Eigen::VectorXd x_vec = b_vector; - - // Fill diagonals vectors - b_diag[0] = A_matrix.coeff(0, 0); - c_diag[0] = A_matrix.coeff(0, 1); - - for (int i = 1; i < n - 1; i++) { - a_diag[i] = A_matrix.coeff(i, i - 1); - b_diag[i] = A_matrix.coeff(i, i); - c_diag[i] = A_matrix.coeff(i, i + 1); - } - a_diag[n - 1] = A_matrix.coeff(n - 1, n - 2); - b_diag[n - 1] = A_matrix.coeff(n - 1, n - 1); - - // start solving - c_diag and x_vec are overwritten - n--; - c_diag[0] /= b_diag[0]; - x_vec[0] /= b_diag[0]; - - for (int i = 1; i < n; i++) { - c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i - 1]; - x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i - 1]) / - (b_diag[i] - a_diag[i] * c_diag[i - 1]); - } - - x_vec[n] = (x_vec[n] - a_diag[n] * x_vec[n - 1]) / - (b_diag[n] - a_diag[n] * c_diag[n - 1]); - - for (int i = n; i-- > 0;) { - x_vec[i] -= c_diag[i] * x_vec[i + 1]; - } - - return x_vec; -} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5f395cd..5210b95 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -11,7 +11,7 @@ if(NOT DOCTEST_LIB) FetchContent_MakeAvailable(DocTest) endif() -add_executable(testTug setup.cpp testSimulation.cpp testGrid.cpp testFTCS.cpp) # testBoundaryCondition.cpp testDiffusion.cpp +add_executable(testTug setup.cpp testSimulation.cpp testGrid.cpp testFTCS.cpp) target_link_libraries(testTug doctest tug) # get relative path of the CSV file diff --git a/test/testBoundaryCondition.cpp b/test/testBoundaryCondition.cpp deleted file mode 100644 index 7b71312..0000000 --- a/test/testBoundaryCondition.cpp +++ /dev/null @@ -1,196 +0,0 @@ -#include -#include - -using namespace tug::bc; - -#define BC_CONST_VALUE 1e-5 - -TEST_CASE("1D Boundary Condition") { - - BoundaryCondition bc(5); - boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE}; - - SUBCASE("valid get") { CHECK_EQ(bc(BC_SIDE_LEFT).value, 0); } - - SUBCASE("invalid get") { - CHECK_THROWS(bc(BC_SIDE_TOP)); - CHECK_THROWS(bc(BC_SIDE_LEFT, 1)); - } - - SUBCASE("valid set") { - CHECK_NOTHROW(bc(BC_SIDE_LEFT) = bc_set); - CHECK_EQ(bc(BC_SIDE_LEFT).value, bc_set.value); - CHECK_EQ(bc(BC_SIDE_LEFT).type, bc_set.type); - } - - SUBCASE("invalid set") { - CHECK_THROWS(bc(BC_SIDE_TOP) = bc_set); - CHECK_THROWS(bc(BC_SIDE_LEFT, 0) = bc_set); - } - - SUBCASE("valid row getter") { - bc(BC_SIDE_LEFT) = bc_set; - bc_tuple tup = bc.row_boundary(0); - - CHECK_EQ(tup[0].value, BC_CONST_VALUE); - CHECK_EQ(tup[1].value, 0); - } - - SUBCASE("invalid row and col getter") { - CHECK_THROWS(bc.row_boundary(1)); - CHECK_THROWS(bc.col_boundary(0)); - } -} - -TEST_CASE("2D Boundary Condition") { - - BoundaryCondition bc(5, 5); - boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE}; - - SUBCASE("valid get") { CHECK_EQ(bc(BC_SIDE_LEFT, 0).value, 0); } - - SUBCASE("invalid get") { - CHECK_THROWS(bc(5, 0)); - CHECK_THROWS(bc(BC_SIDE_LEFT)); - } - - SUBCASE("valid set") { - CHECK_NOTHROW(bc(BC_SIDE_LEFT, 0) = bc_set); - CHECK_EQ(bc(BC_SIDE_LEFT, 0).value, bc_set.value); - CHECK_EQ(bc(BC_SIDE_LEFT, 0).type, bc_set.type); - } - - SUBCASE("invalid set") { - CHECK_THROWS(bc(BC_SIDE_LEFT) = bc_set); - CHECK_THROWS(bc(5, 0) = bc_set); - } - - SUBCASE("call of setSide") { - CHECK_NOTHROW(bc.setSide(BC_SIDE_BOTTOM, bc_set)); - CHECK_EQ(bc(BC_SIDE_BOTTOM, 1).value, bc_set.value); - CHECK_EQ(bc(BC_SIDE_BOTTOM, 1).type, bc_set.type); - } - - SUBCASE("get and set of side") { - std::vector bc_vec; - CHECK_NOTHROW(bc_vec = bc.getSide(BC_SIDE_BOTTOM)); - bc_vec[3] = {BC_TYPE_CONSTANT, 1e-5}; - CHECK_NOTHROW(bc.setSide(BC_SIDE_BOTTOM, bc_vec)); - CHECK_EQ(bc(BC_SIDE_BOTTOM, 3).type, BC_TYPE_CONSTANT); - CHECK_EQ(bc(BC_SIDE_BOTTOM, 3).value, 1e-5); - - CHECK_EQ(bc(BC_SIDE_BOTTOM, 2).value, 0); - } -} - -TEST_CASE("Boundary Condition helpers") { - boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE}; - - SUBCASE("return boundary condition skeleton") { - boundary_condition bc_test = - BoundaryCondition::returnBoundaryCondition(bc_set.type, bc_set.value); - CHECK_EQ(bc_test.value, bc_set.value); - CHECK_EQ(bc_test.type, bc_set.type); - } -} - -TEST_CASE("1D special inner grid cells") { - BoundaryCondition bc(5); - boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE}; - - SUBCASE("valid set") { CHECK_NOTHROW(bc.setInnerBC(bc_set, 0, 0)); } - - SUBCASE("valid get") { - bc.setInnerBC(bc_set, 0, 0); - CHECK_EQ(bc.getInnerBC(0, 0).type, bc_set.type); - } - - SUBCASE("invalid get") { - CHECK_EQ(bc.getInnerBC(1, 0).type, BC_UNSET); - CHECK_THROWS(bc.getInnerBC(0, 1)); - } - - SUBCASE("invalid set") { CHECK_THROWS(bc.setInnerBC(bc_set, 0, 1)); } - - SUBCASE("valid row getter") { - bc.setInnerBC(bc_set, 1, 0); - - bc_vec ret = bc.getInnerRow(0); - - CHECK_EQ(ret[0].type, BC_UNSET); - CHECK_EQ(ret[1].type, bc_set.type); - } - - SUBCASE("invalid row getter") { CHECK_THROWS(bc.getInnerRow(1)); } - - SUBCASE("invalid col getter") { CHECK_THROWS(bc.getInnerCol(0)); } -} - -TEST_CASE("2D special inner grid cells") { - BoundaryCondition bc(5, 5); - boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE}; - - SUBCASE("valid set") { CHECK_NOTHROW(bc.setInnerBC(bc_set, 0, 0)); } - - SUBCASE("valid get") { - bc.setInnerBC(bc_set, 0, 0); - CHECK_EQ(bc.getInnerBC(0, 0).type, bc_set.type); - } - - SUBCASE("invalid get") { - CHECK_EQ(bc.getInnerBC(1, 0).type, BC_UNSET); - CHECK_THROWS(bc.getInnerBC(5, 5)); - } - - SUBCASE("invalid set") { CHECK_THROWS(bc.setInnerBC(bc_set, 5, 5)); } - - SUBCASE("valid row getter") { - bc.setInnerBC(bc_set, 0, 0); - bc.setInnerBC(bc_set, 1, 1); - - bc_vec ret = bc.getInnerRow(0); - - CHECK_EQ(ret[0].type, bc_set.type); - CHECK_EQ(ret[1].type, BC_UNSET); - - ret = bc.getInnerRow(1); - - CHECK_EQ(ret[0].type, BC_UNSET); - CHECK_EQ(ret[1].type, bc_set.type); - } - - SUBCASE("valid col getter") { - bc.setInnerBC(bc_set, 0, 1); - bc.setInnerBC(bc_set, 1, 0); - - bc_vec ret = bc.getInnerCol(0); - - CHECK_EQ(ret[0].type, BC_UNSET); - CHECK_EQ(ret[1].type, bc_set.type); - - ret = bc.getInnerCol(1); - - CHECK_EQ(ret[0].type, bc_set.type); - CHECK_EQ(ret[1].type, BC_UNSET); - } - - SUBCASE("unset boundary condition") { - bc.setInnerBC(bc_set, 0, 0); - bc.setInnerBC(bc_set, 1, 1); - bc.unsetInnerBC(1, 1); - - bc_vec ret = bc.getInnerRow(1); - - CHECK_EQ(ret[0].type, BC_UNSET); - CHECK_EQ(ret[1].type, BC_UNSET); - - ret = bc.getInnerCol(1); - - CHECK_EQ(ret[0].type, BC_UNSET); - CHECK_EQ(ret[1].type, BC_UNSET); - } - - SUBCASE("invalid row getter") { CHECK_THROWS(bc.getInnerRow(5)); } - - SUBCASE("invalid col getter") { CHECK_THROWS(bc.getInnerCol(5)); } -} diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp deleted file mode 100644 index 7039a14..0000000 --- a/test/testDiffusion.cpp +++ /dev/null @@ -1,170 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace tug::bc; -using namespace tug::diffusion; - -#define DIMENSION 2 -#define N 51 -#define M 51 -#define MID 1300 - -static std::vector alpha(N *M, 1e-3); - -static TugInput setupDiffu() { - TugInput diffu; - - diffu.setTimestep(1); - diffu.setGridCellN(N, M); - diffu.setDomainSize(N, M); - - return diffu; -} - -TEST_CASE("closed boundaries - 1 concentration to 1 - rest 0") { - std::vector field(N * M, 0); - - field[MID] = 1; - - BoundaryCondition bc(N, M); - - TugInput diffu = setupDiffu(); - - uint32_t iterations = 1000; - double sum = 0; - - for (int t = 0; t < iterations; t++) { - ADI_2D(diffu, field.data(), alpha.data()); - } - - // iterate through rows - for (int i = 0; i < M; i++) { - // iterate through columns - for (int j = 0; j < N; j++) { - sum += field[i * N + j]; - } - } - CAPTURE(sum); - // epsilon of 1e-8 - CHECK(sum == doctest::Approx(1).epsilon(1e-6)); -} - -TEST_CASE("constant boundaries (0) - 1 concentration to 1 - rest 0") { - std::vector field(N * M, 0); - - field[MID] = 1; - - BoundaryCondition bc(N, M); - - boundary_condition input = {BC_TYPE_CONSTANT, 0}; - - bc.setSide(BC_SIDE_LEFT, input); - bc.setSide(BC_SIDE_RIGHT, input); - bc.setSide(BC_SIDE_TOP, input); - bc.setSide(BC_SIDE_BOTTOM, input); - - TugInput diffu = setupDiffu(); - diffu.setBoundaryCondition(bc); - - uint32_t max_iterations = 20000; - bool reached = false; - - int t = 0; - - for (t = 0; t < max_iterations; t++) { - ADI_2D(diffu, field.data(), alpha.data()); - - if (field[N * M - 1] > 1e-15) { - reached = true; - break; - } - } - - if (!reached) { - CAPTURE(field[N * M - 1]); - FAIL_CHECK( - "Concentration did not reach boundaries after count of iterations: ", - t); - } -} - -TEST_CASE( - "constant top and bottom (1 and 0) - left and right closed - 0 inlet") { - std::vector field(N * M, 0); - - BoundaryCondition bc(N, M); - - boundary_condition top = - BoundaryCondition::returnBoundaryCondition(BC_TYPE_CONSTANT, 1); - boundary_condition bottom = - BoundaryCondition::returnBoundaryCondition(BC_TYPE_CONSTANT, 0); - - bc.setSide(BC_SIDE_TOP, top); - bc.setSide(BC_SIDE_BOTTOM, bottom); - - TugInput diffu = setupDiffu(); - diffu.setBoundaryCondition(bc); - - uint32_t max_iterations = 100; - - for (int t = 0; t < max_iterations; t++) { - ADI_2D(diffu, field.data(), alpha.data()); - } - - for (int i = 0; i < N; i++) { - double above = field[i]; - for (int j = 1; j < M; j++) { - double curr = field[j * N + i]; - if (curr > above) { - CAPTURE(curr); - CAPTURE(above); - FAIL("Concentration below is greater than above @ cell ", j * N + i); - } - } - } -} - -TEST_CASE("2D closed boundaries, 1 constant cell in the middle") { - std::vector field(N * M, 0); - double val = 1e-2; - - BoundaryCondition bc(N, M); - - field[MID] = val; - bc.setInnerBC({BC_TYPE_CONSTANT, val}, N / 2, M / 2); - - TugInput diffu = setupDiffu(); - diffu.setBoundaryCondition(bc); - - uint32_t max_iterations = 100; - - double sum = val; - - for (int t = 0; t < max_iterations; t++) { - ADI_2D(diffu, field.data(), alpha.data()); - - CHECK_EQ(field[MID], val); - - double new_sum = .0; - - for (int i = 0; i < M; i++) { - // iterate through columns - for (int j = 0; j < N; j++) { - new_sum += field[i * N + j]; - } - } - - if (sum > new_sum) { - CAPTURE(sum); - CAPTURE(new_sum); - FAIL("Sum of concentrations is equal or lower than to previous iteration " - "after iteration ", - t + 1); - } - - sum = new_sum; - } -}