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

added Thomas Solver with option to choose solver and cleaned up the repository

Browse Source
This commit is contained in:
philippun 2023-08-23 12:24:35 +02:00
parent 27829a1463
commit 32b05a8a87
33 changed files with 181 additions and 2632 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
examples/CMakeLists.txt
View File

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

59
examples/Rcpp-BTCS-1d.cpp
View File

@ -1,59 +0,0 @@
#include "../include/diffusion/BTCSDiffusion.hpp"
#include "../include/diffusion/BoundaryCondition.hpp"
#include <algorithm> // for copy, max
#include <cmath>
#include <iomanip>
#include <iostream> // for std
#include <vector> // for vector
#include <Rcpp.h>
using namespace std;
using namespace Diffusion;
using namespace Rcpp;
//using namespace Eigen;
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
std::vector<double> & diff1D(int n,
double length,
std::vector<double> & field,
std::vector<double> & 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<double> alpha(n, 1 * pow(10, -1));
// std::vector<double> field(n, 1 * std::pow(10, -6));
std::vector<boundary_condition> bc(n, {0,0});
// create instance of diffusion module
BTCSDiffusion diffu(dim);
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);
}

60
examples/Rcpp-BTCS-2d.cpp
View File

@ -1,60 +0,0 @@
#include "../include/diffusion/BTCSDiffusion.hpp"
#include "../include/diffusion/BoundaryCondition.hpp"
#include <algorithm> // for copy, max
#include <cmath>
#include <iomanip>
#include <iostream> // for std
#include <vector> // for vector
#include <Rcpp.h>
using namespace std;
using namespace Diffusion;
using namespace Rcpp;
//using namespace Eigen;
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
std::vector<double> & diff2D(int nx,
int ny,
double lenx,
double leny,
std::vector<double> & field,
std::vector<double> & alpha,
double timestep,
int iterations)
{
// problem dimensionality
int dim = 2;
// total number of grid cells
int n = nx*ny;
std::vector<boundary_condition> 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);
}

159
examples/Rcpp-interface.R
View File

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

59
examples/Rcpp-interface.cpp
View File

@ -1,59 +0,0 @@
#include "../include/diffusion/BTCSDiffusion.hpp"
#include "../include/diffusion/BoundaryCondition.hpp"
#include <algorithm> // for copy, max
#include <cmath>
#include <iomanip>
#include <iostream> // for std
#include <vector> // for vector
#include <Rcpp.h>
using namespace std;
using namespace Diffusion;
using namespace Rcpp;
//using namespace Eigen;
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
std::vector<double> & diff1D(int n,
double length,
std::vector<double> & field,
std::vector<double> & 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<double> alpha(n, 1 * pow(10, -1));
// std::vector<double> field(n, 1 * std::pow(10, -6));
std::vector<boundary_condition> bc(n, {0,0});
// create instance of diffusion module
BTCSDiffusion diffu(dim);
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);
}

40
examples/RcppFTCS.cpp
View File

@ -1,40 +0,0 @@
// Time-stamp: "Last modified 2022-03-15 18:15:39 delucia"
#include <Rcpp.h>
#include <iostream> // for std
#include <vector> // 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);
}

35
examples/boundary_example1D.cpp
View File

@ -1,35 +0,0 @@
#include <algorithm>
#include <tug/BoundaryCondition.hpp>
#include <tug/Diffusion.hpp>
#include <vector>
#include <iostream>
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<boundary_condition> result_left = example.getSide(side);
// vector<boundary_condition> 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;
}

72
examples/first_example.cpp
View File

@ -1,72 +0,0 @@
#include "tug/BoundaryCondition.hpp"
#include <tug/Diffusion.hpp>
#include <iostream>
#include <fstream>
#include <string>
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<double> alpha(n, 1e-1);
// double alpha = 1e-1;
vector<double> 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();
}

55
examples/main_1D.cpp
View File

@ -1,55 +0,0 @@
#include "diffusion/BTCSBoundaryCondition.hpp"
#include <diffusion/BTCSDiffusion.hpp>
#include <iomanip>
#include <iostream> // for std
#include <vector> // for vector
using namespace std;
using namespace Diffusion;
int main(int argc, char *argv[]) {
// dimension of grid
int dim = 1;
int n = 20;
// create input + diffusion coefficients for each grid cell
std::vector<double> alpha(n, 1e-1);
std::vector<double> 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;
}

58
examples/main_2D.cpp
View File

@ -1,58 +0,0 @@
#include "diffusion/BTCSBoundaryCondition.hpp"
#include <diffusion/BTCSDiffusion.hpp>
#include <iomanip>
#include <iostream> // for std
#include <vector> // for vector
using namespace std;
using namespace Diffusion;
int main(int argc, char *argv[]) {
// dimension of grid
int dim = 2;
int n = 5;
int m = 5;
// create input + diffusion coefficients for each grid cell
std::vector<double> alpha(n * m, 1e-6);
std::vector<double> 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;
}

59
examples/main_2D_const.cpp
View File

@ -1,59 +0,0 @@
#include "diffusion/BTCSBoundaryCondition.hpp"
#include <diffusion/BTCSDiffusion.hpp>
#include <iomanip>
#include <iostream> // for std
#include <vector> // for vector
using namespace std;
using namespace Diffusion;
int main(int argc, char *argv[]) {
// dimension of grid
int dim = 2;
int n = 5;
int m = 5;
// create input + diffusion coefficients for each grid cell
std::vector<double> alpha(n * m, 1e-1);
std::vector<double> 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;
}

67
examples/main_2D_mdl.cpp
View File

@ -1,67 +0,0 @@
#include "diffusion/BTCSBoundaryCondition.hpp"
#include <diffusion/BTCSDiffusion.hpp>
#include <iomanip>
#include <iostream> // for std
#include <vector> // for vector
using namespace std;
using namespace Diffusion;
int main(int argc, char *argv[]) {
// dimension of grid
int dim = 2;
int n = 501;
int m = 501;
// create input + diffusion coefficients for each grid cell
std::vector<double> alpha(n * m, 1e-3);
std::vector<double> 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;
}

38
examples/meeting_example.cpp
View File

@ -1,38 +0,0 @@
#include "tug/Boundary.hpp"
#include <tug/Simulation.hpp>
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();
}

50
examples/profiling.cpp
View File

@ -1,50 +0,0 @@
#include <tug/Simulation.hpp>
#include <iostream>
#include <fstream>
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<std::chrono::milliseconds>(end - begin);
myfile << milliseconds.count() << endl;
}
}
cout << endl;
myfile << endl;
}
myfile.close();
}

5
examples/profiling_openmp.cpp
View File

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

25
examples/pseudo_example.cpp
View File

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

93
examples/second_example.cpp
View File

@ -1,93 +0,0 @@
#include "tug/BoundaryCondition.hpp"
#include <tug/Diffusion.hpp>
#include <iostream>
#include <fstream>
#include <string>
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<double> alpha(n * m, 1);
vector<double> 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<n; i++) {
for (int j = 0; j<m; j++) {
cout << field[n * i + j] << " ";
}
cout << endl;
}
cout << endl;
TugInput input_param;
input_param.setTimestep(0.1);
input_param.setGridCellN(n, m);
input_param.setDomainSize(n, m);
BoundaryCondition bc(n, m);
boundary_condition bc_constant;
bc_constant.type = BC_TYPE_CONSTANT;
bc_constant.value = 0;
bc.setSide(BC_SIDE_LEFT, bc_constant);
bc.setSide(BC_SIDE_RIGHT, bc_constant);
bc.setSide(BC_SIDE_BOTTOM, bc_constant);
// bc_constant.value = 2000;
bc.setSide(BC_SIDE_TOP, bc_constant);
input_param.setBoundaryCondition(bc);
// int iterations = 1000;
// for (int t = 0; t < iterations; t++) {
// double result = ADI_2D(input_param, &field[0], &alpha[0]);
// if (t % 100 == 0) {
// cout << "Iteration " << t << ":" << endl;
// for (int i = 0; i<n; i++) {
// for (int j = 0; j<m; j++) {
// cout << field[n * i + j] << " ";
// }
// cout << endl;
// }
// cout << endl;
// }
// }
ofstream myfile;
myfile.open("output.csv");
if (!myfile) {
exit(1);
}
int iterations = 100;
for (int t = 0; t < iterations; t++) {
double result = ADI_2D(input_param, &field[0], &alpha[0]);
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) {
myfile << field[i * m + j];
if (j < m-1) {
myfile << ", ";
}
}
myfile << "\n";
}
myfile << std::endl << std::endl;
}
myfile.close();
}

283
include/tug/BoundaryCondition.hpp
View File

@ -1,283 +0,0 @@
#ifndef BOUNDARYCONDITION_H_
#define BOUNDARYCONDITION_H_
#include <array>
#include <map>
#include <stdexcept>
#include <stdint.h>
#include <vector>
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<boundary_condition, 2> bc_tuple;
typedef std::vector<boundary_condition> 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<boundary_condition> &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<boundary_condition>;
/**
* 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<boundary_condition> bc_internal;
std::map<uint32_t, boundary_condition> inner_cells;
uint8_t dim;
std::array<uint32_t, 2> 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_

140
include/tug/Diffusion.hpp
View File

@ -1,140 +0,0 @@
#ifndef DIFFUSION_H_
#define DIFFUSION_H_
#include "BoundaryCondition.hpp"
#include "Solver.hpp"
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <array>
#include <vector>
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<uint32_t, MAX_ARR_SIZE>
grid_cells; /**< Count of grid cells in each of the 3 directions.*/
std::array<double, MAX_ARR_SIZE>
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<double> &,
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<double> &,
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_

21
include/tug/Simulation.hpp
View File

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

41
include/tug/Solver.hpp
View File

@ -1,41 +0,0 @@
#ifndef SOLVER_H_
#define SOLVER_H_
#include <Eigen/Dense>
#include <Eigen/Sparse>
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<double> &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<double> &A_matrix,
const Eigen::VectorXd &b_vector) -> Eigen::VectorXd;
} // namespace solver
} // namespace tug
#endif // SOLVER_H_

191
include/tug/progressbar.hpp
View File

@ -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 <iostream>
#include <ostream>
#include <string>
#include <stdexcept>
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

328
src/BTCS.cpp
View File

@ -1,328 +0,0 @@
#include <array>
#include <iostream>
#include <tug/BoundaryCondition.hpp>
#include <tug/Diffusion.hpp>
#include <tug/Solver.hpp>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/src/Core/Matrix.h>
#include <chrono>
#include <vector>
#include "TugUtils.hpp"
#ifdef _OPENMP
#include <omp.h>
#else
#define omp_get_thread_num() 0
#endif
inline auto
init_delta(const std::array<double, tug::diffusion::MAX_ARR_SIZE> &domain_size,
const std::array<uint32_t, tug::diffusion::MAX_ARR_SIZE> &grid_cells,
const uint8_t dim) -> std::vector<double> {
std::vector<double> 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<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
using DVectorRowMajor =
Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>;
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<double, 3> 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<double> {
Eigen::SparseMatrix<double> 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<double> &,
const Eigen::VectorXd &)) -> DVectorRowMajor {
const Eigen::SparseMatrix<double> 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<DVectorRowMajor> c_in(field, size);
Eigen::Map<const DVectorRowMajor> 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<DMatrixRowMajor> c_in(field, n_rows, n_cols);
Eigen::Map<const DMatrixRowMajor> 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);
}

76
src/BTCSv2.cpp
View File

@ -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<double> &A, VectorXd &b) {
// use of EigenLU solver
static VectorXd EigenLUAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
SparseLU<SparseMatrix<double>> solver;
solver.analyzePattern(A);
@ -273,9 +274,53 @@ static VectorXd solve(SparseMatrix<double> &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<double> &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 &timestep) {
static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solverFunc) (SparseMatrix<double> &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 &timestep) {
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 &timestep) {
// BTCS solution for 2D grid
static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solverFunc) (SparseMatrix<double> &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 &timestep) {
SparseLU<SparseMatrix<double>> 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 &timestep) {
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 &timestep) {
}
// entry point; differentiate between 1D and 2D grid
static void BTCS(Grid &grid, Boundary &bc, double &timestep) {
// entry point for EigenLU solver; differentiate between 1D and 2D grid
static void BTCS_LU(Grid &grid, Boundary &bc, double &timestep) {
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 &timestep) {
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!");
}

3
src/Boundary.cpp
View File

@ -1,5 +1,4 @@
#include "TugUtils.hpp"
#include "tug/BoundaryCondition.hpp"
#include <iostream>
#include <omp.h>
#include <tug/Boundary.hpp>
@ -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;
}

209
src/BoundaryCondition.cpp
View File

@ -1,209 +0,0 @@
#include <algorithm>
#include <tug/BoundaryCondition.hpp>
#include <vector>
#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<tug::bc::boundary_condition> &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<tug::bc::boundary_condition> {
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<tug::bc::boundary_condition> 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};
}

2
src/CMakeLists.txt
View File

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

9
src/README.md
View File

@ -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
</pre>
@ -23,10 +20,6 @@ src/
**Boundary.cpp** Implementation of Boundary class, that holds the boundary conditions.
**BoundaryCondition.cpp** <i>Outdated</i> implementation of boundary conditions.
**BTCS.cpp** <i>Outdated</i> 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** <i>Outdated</i> implementation of Eigen solvers.
**TugUtils.hpp** Helper functions for other source files.

132
src/Simulation.cpp
View File

@ -7,8 +7,10 @@
#include <fstream>
#ifndef SIMULATION_H_
#define SIMULATION_H_
#include "BTCSv2.cpp"
#include <tug/progressbar.hpp>
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

61
src/Solver.cpp
View File

@ -1,61 +0,0 @@
#include <tug/Solver.hpp>
#include <Eigen/Dense>
#include <Eigen/Sparse>
auto tug::solver::EigenLU(const Eigen::SparseMatrix<double> &A_matrix,
const Eigen::VectorXd &b_vector) -> Eigen::VectorXd {
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
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<double> &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;
}

2
test/CMakeLists.txt
View File

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

196
test/testBoundaryCondition.cpp
View File

@ -1,196 +0,0 @@
#include <doctest/doctest.h>
#include <tug/BoundaryCondition.hpp>
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<boundary_condition> 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)); }
}

170
test/testDiffusion.cpp
View File

@ -1,170 +0,0 @@
#include <doctest/doctest.h>
#include <tug/BoundaryCondition.hpp>
#include <tug/Diffusion.hpp>
#include <tug/Solver.hpp>
#include <vector>
using namespace tug::bc;
using namespace tug::diffusion;
#define DIMENSION 2
#define N 51
#define M 51
#define MID 1300
static std::vector<double> 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<double> 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<double> 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<double> 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<double> 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;
}
}