mirror of
https://git.gfz-potsdam.de/naaice/tug.git
synced 2025-12-13 09:28:23 +01:00
Merge branch 'hannes-philipp' of git.gfz-potsdam.de:naaice/tug into hannes-philipp
This commit is contained in:
Hannes Signer
commit
dde7fb4783
@ -35,10 +35,14 @@ int main(int argc, char *argv[]) {
|
||||
|
||||
// create a boundary with constant values
|
||||
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.setBoundarySideClosed(BC_SIDE_LEFT);
|
||||
bc.setBoundarySideClosed(BC_SIDE_RIGHT);
|
||||
bc.setBoundarySideClosed(BC_SIDE_TOP);
|
||||
bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
|
||||
// bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
|
||||
// bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
|
||||
// bc.setBoundarySideConstant(BC_SIDE_TOP, 0);
|
||||
// bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
|
||||
|
||||
|
||||
// (optional) set boundary condition values for one side, e.g.:
|
||||
@ -66,7 +70,7 @@ int main(int argc, char *argv[]) {
|
||||
simulation.setIterations(100);
|
||||
|
||||
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
|
||||
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
|
||||
simulation.setOutputCSV(CSV_OUTPUT_XTREME);
|
||||
|
||||
// **** RUN SIMULATION ****
|
||||
|
||||
|
||||
@ -78,6 +78,7 @@ int main(int argc, char *argv[]) {
|
||||
|
||||
// set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
|
||||
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
|
||||
|
||||
|
||||
// **** RUN SIMULATION ****
|
||||
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
/**
|
||||
* @file Boundary.hpp
|
||||
* @brief
|
||||
*
|
||||
* @brief API of Boundary class, that holds all information for each boundary condition
|
||||
* at the edges of the diffusion grid.
|
||||
*
|
||||
*/
|
||||
#ifndef BOUNDARY_H_
|
||||
@ -13,11 +13,19 @@
|
||||
using namespace std;
|
||||
using namespace Eigen;
|
||||
|
||||
/**
|
||||
* @brief Enum defining the two implemented boundary conditions.
|
||||
*
|
||||
*/
|
||||
enum BC_TYPE {
|
||||
BC_TYPE_CLOSED,
|
||||
BC_TYPE_CONSTANT
|
||||
};
|
||||
|
||||
/**
|
||||
* @brief Enum defining all 4 possible sides to a 1D and 2D grid.
|
||||
*
|
||||
*/
|
||||
enum BC_SIDE {
|
||||
BC_SIDE_LEFT,
|
||||
BC_SIDE_RIGHT,
|
||||
@ -149,12 +157,18 @@ class Boundary {
|
||||
* @brief Returns the boundary condition of a specified side as a vector
|
||||
* of BoundarsElement objects.
|
||||
*
|
||||
* @param side Boundary side from which the boundaryconditions are to be returned.
|
||||
* @param side Boundary side from which the boundary conditions are to be returned.
|
||||
* @return vector<BoundaryElement> Contains the boundary conditions as BoundaryElement objects.
|
||||
*/
|
||||
vector<BoundaryElement> getBoundarySide(BC_SIDE side);
|
||||
|
||||
// TODO write documentation and tests for this method
|
||||
/**
|
||||
* @brief Get thes Boundary Side Values as a vector. Value is -1 in case some specific
|
||||
boundary is closed.
|
||||
*
|
||||
* @param side Boundary side for which the values are to be returned.
|
||||
* @return VectorXd Vector with values as doubles.
|
||||
*/
|
||||
VectorXd getBoundarySideValues(BC_SIDE side);
|
||||
|
||||
/**
|
||||
@ -192,9 +206,9 @@ class Boundary {
|
||||
double getBoundaryElementValue(BC_SIDE side, int index);
|
||||
|
||||
private:
|
||||
Grid grid;
|
||||
Grid grid; // Boundary is directly dependent on the dimensions of a predefined
|
||||
|
||||
vector<vector<BoundaryElement>> boundaries;
|
||||
vector<vector<BoundaryElement>> boundaries; // Vector with Boundary Element information
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
@ -132,6 +132,13 @@ class Grid {
|
||||
*/
|
||||
void setDomain(int domainRow, int domainCol);
|
||||
|
||||
/**
|
||||
* @brief Gets the delta value for 1D-Grid. Grid must be one dimensional.
|
||||
*
|
||||
* @return double Delta value.
|
||||
*/
|
||||
double getDelta();
|
||||
|
||||
/**
|
||||
* @brief Gets the delta value in x-direction.
|
||||
*
|
||||
@ -156,8 +163,8 @@ class Grid {
|
||||
int domainRow; // number of domain rows
|
||||
double deltaCol; // delta in x-direction (between columns)
|
||||
double deltaRow; // delta in y-direction (between rows)
|
||||
MatrixXd concentrations;
|
||||
MatrixXd alphaX;
|
||||
MatrixXd alphaY;
|
||||
MatrixXd concentrations; // Matrix holding grid concentrations
|
||||
MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction
|
||||
MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction
|
||||
|
||||
};
|
||||
@ -1,34 +1,52 @@
|
||||
/**
|
||||
* @file Simulation.hpp
|
||||
* @brief
|
||||
* @brief API of Simulation class, that holds all information regarding a specific simulation
|
||||
* run like its timestep, number of iterations and output options. Simulation object
|
||||
* also holds a predefined Grid and Boundary object.
|
||||
*
|
||||
*/
|
||||
#include "Boundary.hpp"
|
||||
#include <ios>
|
||||
|
||||
using namespace std;
|
||||
|
||||
/**
|
||||
* @brief Enum defining the two implemented solution approaches.
|
||||
*
|
||||
*/
|
||||
enum APPROACH {
|
||||
FTCS_APPROACH, // Forward Time-Centered Space
|
||||
BTCS_APPROACH // Backward Time-Centered Space
|
||||
};
|
||||
|
||||
/**
|
||||
* @brief Enum holding different options for .csv output.
|
||||
*
|
||||
*/
|
||||
enum CSV_OUTPUT {
|
||||
CSV_OUTPUT_OFF, // do not produce csv output
|
||||
CSV_OUTPUT_ON, // produce csv output with last concentration matrix
|
||||
CSV_OUTPUT_VERBOSE, // produce csv output with all concentration matrices
|
||||
CSV_OUTPUT_XTREME // produce csv output with all concentration matrices and simulation environment
|
||||
CSV_OUTPUT_XTREME // csv output like VERBOSE but additional boundary conditions at beginning
|
||||
};
|
||||
|
||||
/**
|
||||
* @brief Enum holding different options for console output.
|
||||
*
|
||||
*/
|
||||
enum CONSOLE_OUTPUT {
|
||||
CONSOLE_OUTPUT_OFF, // do not print any output to console
|
||||
CONSOLE_OUTPUT_ON, // print before and after concentrations to console
|
||||
CONSOLE_OUTPUT_VERBOSE // print all concentration matrices to console
|
||||
};
|
||||
|
||||
/**
|
||||
* @brief Enum holding different options for time measurement.
|
||||
*
|
||||
*/
|
||||
enum TIME_MEASURE {
|
||||
TIME_MEASURE_OFF, // do not print any time measures
|
||||
TIME_MEASURE_ON, // print time measure after last iteration
|
||||
TIME_MEASURE_VERBOSE // print time measures after each iteration
|
||||
TIME_MEASURE_ON // print time measure after last iteration
|
||||
};
|
||||
|
||||
/**
|
||||
|
||||
256
src/BTCSv2.cpp
256
src/BTCSv2.cpp
@ -1,19 +1,86 @@
|
||||
/**
|
||||
* @file BTCS.cpp
|
||||
* @brief Implementation of heterogenous BTCS (backward time-centered space) solution
|
||||
* of diffusion equation in 1D and 2D space.
|
||||
*
|
||||
*/
|
||||
|
||||
#include "FTCS.cpp"
|
||||
#include <tug/Boundary.hpp>
|
||||
|
||||
using namespace Eigen;
|
||||
|
||||
|
||||
static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, int numCols, int rowIndex, double sx) {
|
||||
// calculates coefficient for left boundary in constant case
|
||||
static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, double &sx) {
|
||||
double centerCoeff;
|
||||
double rightCoeff;
|
||||
|
||||
centerCoeff = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1))
|
||||
+ 2 * alpha(rowIndex,0));
|
||||
rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1));
|
||||
|
||||
return {centerCoeff, rightCoeff};
|
||||
}
|
||||
|
||||
|
||||
// calculates coefficient for left boundary in closed case
|
||||
static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, double &sx) {
|
||||
double centerCoeff;
|
||||
double rightCoeff;
|
||||
|
||||
centerCoeff = 1 + sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1));
|
||||
rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1));
|
||||
|
||||
return {centerCoeff, rightCoeff};
|
||||
}
|
||||
|
||||
|
||||
// calculates coefficient for right boundary in constant case
|
||||
static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, int &n, double &sx) {
|
||||
double leftCoeff;
|
||||
double centerCoeff;
|
||||
|
||||
leftCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n));
|
||||
centerCoeff = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n))
|
||||
+ 2 * alpha(rowIndex,n));
|
||||
|
||||
return {leftCoeff, centerCoeff};
|
||||
}
|
||||
|
||||
|
||||
// calculates coefficient for right boundary in closed case
|
||||
static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, int &n, double &sx) {
|
||||
double leftCoeff;
|
||||
double centerCoeff;
|
||||
|
||||
leftCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n));
|
||||
centerCoeff = 1 + sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n));
|
||||
|
||||
return {leftCoeff, centerCoeff};
|
||||
}
|
||||
|
||||
|
||||
// creates coefficient matrix for next time step from alphas in x-direction
|
||||
static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryElement> bcLeft, vector<BoundaryElement> bcRight, int numCols, int rowIndex, double sx) {
|
||||
|
||||
// square matrix of column^2 dimension for the coefficients
|
||||
SparseMatrix<double> cm(numCols, numCols);
|
||||
cm.reserve(VectorXi::Constant(numCols, 3));
|
||||
|
||||
// left column
|
||||
cm.insert(0,0) = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1))
|
||||
+ 2 * alpha(rowIndex,0));
|
||||
cm.insert(0,1) = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1));
|
||||
BC_TYPE type = bcLeft[rowIndex].getType();
|
||||
if (type == BC_TYPE_CONSTANT) {
|
||||
auto [centerCoeffTop, rightCoeffTop] = calcLeftBoundaryCoeffConstant(alpha, rowIndex, sx);
|
||||
cm.insert(0,0) = centerCoeffTop;
|
||||
cm.insert(0,1) = rightCoeffTop;
|
||||
} else if (type == BC_TYPE_CLOSED) {
|
||||
auto [centerCoeffTop, rightCoeffTop] = calcLeftBoundaryCoeffClosed(alpha, rowIndex, sx);
|
||||
cm.insert(0,0) = centerCoeffTop;
|
||||
cm.insert(0,1) = rightCoeffTop;
|
||||
} else {
|
||||
throw_invalid_argument("Undefined Boundary Condition Type somewhere on Left or Top!");
|
||||
}
|
||||
|
||||
// inner columns
|
||||
int n = numCols-1;
|
||||
@ -28,22 +95,106 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, int numCols, int
|
||||
}
|
||||
|
||||
// right column
|
||||
cm.insert(n,n-1) = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n));
|
||||
cm.insert(n,n) = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n))
|
||||
+ 2 * alpha(rowIndex,n));
|
||||
type = bcRight[rowIndex].getType();
|
||||
if (type == BC_TYPE_CONSTANT) {
|
||||
auto [leftCoeffBottom, centerCoeffBottom] = calcRightBoundaryCoeffConstant(alpha, rowIndex, n, sx);
|
||||
cm.insert(n,n-1) = leftCoeffBottom;
|
||||
cm.insert(n,n) = centerCoeffBottom;
|
||||
} else if (type == BC_TYPE_CLOSED) {
|
||||
auto [leftCoeffBottom, centerCoeffBottom] = calcRightBoundaryCoeffClosed(alpha, rowIndex, n, sx);
|
||||
cm.insert(n,n-1) = leftCoeffBottom;
|
||||
cm.insert(n,n) = centerCoeffBottom;
|
||||
} else {
|
||||
throw_invalid_argument("Undefined Boundary Condition Type somewhere on Right or Bottom!");
|
||||
}
|
||||
|
||||
cm.makeCompressed();
|
||||
cm.makeCompressed(); // important for Eigen solver
|
||||
|
||||
return cm;
|
||||
}
|
||||
|
||||
|
||||
// calculates explicity concentration at top boundary in constant case
|
||||
static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrations,
|
||||
MatrixXd &alpha, vector<BoundaryElement> &bcTop, int &rowIndex, int &i, double &sy) {
|
||||
double c;
|
||||
|
||||
c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
|
||||
* concentrations(rowIndex,i)
|
||||
+ (
|
||||
1 - sy * (
|
||||
calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
|
||||
+ 2 * alpha(rowIndex,i)
|
||||
)
|
||||
) * concentrations(rowIndex,i)
|
||||
+ sy * alpha(rowIndex,i) * bcTop[i].getValue();
|
||||
|
||||
return c;
|
||||
}
|
||||
|
||||
|
||||
// calculates explicit concentration at top boundary in closed case
|
||||
static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentrations,
|
||||
MatrixXd &alpha, int &rowIndex, int &i, double &sy) {
|
||||
double c;
|
||||
|
||||
c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
|
||||
* concentrations(rowIndex,i)
|
||||
+ (
|
||||
1 - sy * (
|
||||
calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
|
||||
)
|
||||
) * concentrations(rowIndex,i);
|
||||
|
||||
return c;
|
||||
}
|
||||
|
||||
|
||||
// calculates explicit concentration at bottom boundary in constant case
|
||||
static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concentrations,
|
||||
MatrixXd &alpha, vector<BoundaryElement> &bcBottom, int &rowIndex, int &i, double &sy) {
|
||||
double c;
|
||||
|
||||
c = sy * alpha(rowIndex,i) * bcBottom[i].getValue()
|
||||
+ (
|
||||
1 - sy * (
|
||||
2 * alpha(rowIndex,i)
|
||||
+ calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i))
|
||||
)
|
||||
) * concentrations(rowIndex,i)
|
||||
+ sy * calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i))
|
||||
* concentrations(rowIndex-1,i);
|
||||
|
||||
return c;
|
||||
}
|
||||
|
||||
|
||||
// calculates explicit concentration at bottom boundary in closed case
|
||||
static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentrations,
|
||||
MatrixXd &alpha, int &rowIndex, int &i, double &sy) {
|
||||
double c;
|
||||
|
||||
c = (
|
||||
1 - sy * (
|
||||
+ calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i))
|
||||
)
|
||||
) * concentrations(rowIndex,i)
|
||||
+ sy * calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i))
|
||||
* concentrations(rowIndex-1,i);
|
||||
|
||||
return c;
|
||||
}
|
||||
|
||||
|
||||
// creates a solution vector for next time step from the current state of concentrations
|
||||
static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY,
|
||||
VectorXd &bcLeft, VectorXd &bcRight, VectorXd &bcTop,
|
||||
VectorXd &bcBottom, int length, int rowIndex, double sx, double sy) {
|
||||
vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight,
|
||||
vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom,
|
||||
int length, int rowIndex, double sx, double sy) {
|
||||
|
||||
VectorXd sv(length);
|
||||
int numRows = concentrations.rows();
|
||||
BC_TYPE type;
|
||||
|
||||
// inner rows
|
||||
if (rowIndex > 0 && rowIndex < numRows-1) {
|
||||
@ -65,45 +216,47 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
|
||||
// first row
|
||||
if (rowIndex == 0) {
|
||||
for (int i = 0; i < length; i++) {
|
||||
sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i))
|
||||
* concentrations(rowIndex,i)
|
||||
+ (
|
||||
1 - sy * (
|
||||
calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i))
|
||||
+ 2 * alphaY(rowIndex,i)
|
||||
)
|
||||
) * concentrations(rowIndex,i)
|
||||
+ sy * alphaY(rowIndex,i) * bcTop(i)
|
||||
;
|
||||
type = bcTop[i].getType();
|
||||
if (type == BC_TYPE_CONSTANT) {
|
||||
sv(i) = calcExplicitConcentrationsTopBoundaryConstant(concentrations, alphaY, bcTop, rowIndex, i, sy);
|
||||
} else if (type == BC_TYPE_CLOSED) {
|
||||
sv(i) = calcExplicitConcentrationsTopBoundaryClosed(concentrations, alphaY, rowIndex, i, sy);
|
||||
} else {
|
||||
throw_invalid_argument("Undefined Boundary Condition Type somewhere on Left or Top!");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// last row
|
||||
if (rowIndex == numRows-1) {
|
||||
for (int i = 0; i < length; i++) {
|
||||
sv(i) = sy * alphaY(rowIndex,i) * bcBottom(i)
|
||||
+ (
|
||||
1 - sy * (
|
||||
2 * alphaY(rowIndex,i)
|
||||
+ calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i))
|
||||
)
|
||||
) * concentrations(rowIndex,i)
|
||||
+ sy * calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i))
|
||||
* concentrations(rowIndex-1,i)
|
||||
;
|
||||
type = bcBottom[i].getType();
|
||||
if (type == BC_TYPE_CONSTANT) {
|
||||
sv(i) = calcExplicitConcentrationsBottomBoundaryConstant(concentrations, alphaY, bcBottom, rowIndex, i, sy);
|
||||
} else if (type == BC_TYPE_CLOSED) {
|
||||
sv(i) = calcExplicitConcentrationsBottomBoundaryClosed(concentrations, alphaY, rowIndex, i, sy);
|
||||
} else {
|
||||
throw_invalid_argument("Undefined Boundary Condition Type somewhere on Right or Bottom!");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// first column -> additional fixed concentration change from perpendicular dimension
|
||||
sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft(rowIndex);
|
||||
// first column -> additional fixed concentration change from perpendicular dimension in constant bc case
|
||||
if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) {
|
||||
sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft[rowIndex].getValue();
|
||||
}
|
||||
|
||||
// last column -> additional fixed concentration change from perpendicular dimension
|
||||
sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight(rowIndex);
|
||||
// last column -> additional fixed concentration change from perpendicular dimension in constant bc case
|
||||
if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) {
|
||||
sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight[rowIndex].getValue();
|
||||
}
|
||||
|
||||
return sv;
|
||||
}
|
||||
|
||||
|
||||
// solver for linear equation system; A corresponds to coefficient matrix,
|
||||
// b to the solution vector
|
||||
static VectorXd solve(SparseMatrix<double> &A, VectorXd &b) {
|
||||
SparseLU<SparseMatrix<double>> solver;
|
||||
solver.analyzePattern(A);
|
||||
@ -116,7 +269,7 @@ static VectorXd solve(SparseMatrix<double> &A, VectorXd &b) {
|
||||
// BTCS solution for 1D grid
|
||||
static void BTCS_1D(Grid &grid, Boundary &bc, double ×tep) {
|
||||
int length = grid.getLength();
|
||||
double sx = timestep / (grid.getDeltaCol() * grid.getDeltaCol()); // TODO create method getDelta() for 1D case
|
||||
double sx = timestep / (grid.getDelta() * grid.getDelta());
|
||||
|
||||
VectorXd concentrations_t1(length);
|
||||
|
||||
@ -124,16 +277,20 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double ×tep) {
|
||||
VectorXd b(length);
|
||||
|
||||
MatrixXd alpha = grid.getAlpha();
|
||||
VectorXd bcLeft = bc.getBoundarySideValues(BC_SIDE_LEFT);
|
||||
VectorXd bcRight = bc.getBoundarySideValues(BC_SIDE_RIGHT);
|
||||
vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
|
||||
vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
|
||||
|
||||
MatrixXd concentrations = grid.getConcentrations();
|
||||
A = createCoeffMatrix(alpha, length, 0, sx); // this is exactly same as in 2D
|
||||
A = createCoeffMatrix(alpha, bcLeft, bcRight, length, 0, sx); // this is exactly same as in 2D
|
||||
for (int i = 0; i < length; i++) {
|
||||
b(i) = concentrations(0,i); // TODO check if this is really the only thing on right hand side of equation in 1D?
|
||||
b(i) = concentrations(0,i);
|
||||
}
|
||||
if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) {
|
||||
b(0) += 2 * sx * alpha(0,0) * bcLeft[0].getValue();
|
||||
}
|
||||
if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) {
|
||||
b(length-1) += 2 * sx * alpha(0,length-1) * bcRight[0].getValue();
|
||||
}
|
||||
b(0) += 2 * sx * alpha(0,0) * bcLeft(0);
|
||||
b(length-1) += 2 * sx * alpha(0,length-1) * bcRight(0);
|
||||
|
||||
concentrations_t1 = solve(A, b);
|
||||
|
||||
@ -160,21 +317,21 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep) {
|
||||
|
||||
MatrixXd alphaX = grid.getAlphaX();
|
||||
MatrixXd alphaY = grid.getAlphaY();
|
||||
VectorXd bcLeft = bc.getBoundarySideValues(BC_SIDE_LEFT);
|
||||
VectorXd bcRight = bc.getBoundarySideValues(BC_SIDE_RIGHT);
|
||||
VectorXd bcTop = bc.getBoundarySideValues(BC_SIDE_TOP);
|
||||
VectorXd bcBottom = bc.getBoundarySideValues(BC_SIDE_BOTTOM);
|
||||
vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
|
||||
vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
|
||||
vector<BoundaryElement> bcTop = bc.getBoundarySide(BC_SIDE_TOP);
|
||||
vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
|
||||
|
||||
MatrixXd concentrations = grid.getConcentrations();
|
||||
for (int i = 0; i < rowMax; i++) {
|
||||
|
||||
A = createCoeffMatrix(alphaX, colMax, i, sx);
|
||||
A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx);
|
||||
b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
|
||||
bcTop, bcBottom, colMax, i, sx, sy);
|
||||
row_t1 = solve(A, b);
|
||||
|
||||
for (int j = 0; j < colMax; j++) {
|
||||
concentrations_t1(i,j) = row_t1(j);
|
||||
concentrations_t1(i,j) = row_t1(j); // can potentially be improved by using Eigen method
|
||||
}
|
||||
|
||||
}
|
||||
@ -184,13 +341,14 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep) {
|
||||
alphaY.transposeInPlace();
|
||||
for (int i = 0; i < colMax; i++) {
|
||||
|
||||
A = createCoeffMatrix(alphaY, rowMax, i, sy);
|
||||
// swap alphas, boundary conditions and sx/sy for column-wise calculation
|
||||
A = createCoeffMatrix(alphaY, bcLeft, bcRight, rowMax, i, sy);
|
||||
b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
|
||||
bcLeft, bcRight, rowMax, i, sy, sx);
|
||||
row_t1 = solve(A, b);
|
||||
|
||||
for (int j = 0; j < rowMax; j++) {
|
||||
concentrations(i,j) = row_t1(j);
|
||||
concentrations(i,j) = row_t1(j); // can potentially be improved by using Eigen method
|
||||
}
|
||||
}
|
||||
concentrations.transposeInPlace();
|
||||
|
||||
@ -115,6 +115,10 @@ 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) {
|
||||
values(i) = -1;
|
||||
continue;
|
||||
}
|
||||
values(i) = getBoundaryElementValue(side, i);
|
||||
}
|
||||
|
||||
|
||||
10
src/Grid.cpp
10
src/Grid.cpp
@ -13,7 +13,6 @@ Grid::Grid(int length) {
|
||||
this->deltaCol = double(this->domainCol)/double(this->col); // -> 1
|
||||
this->dim = 1;
|
||||
|
||||
// TODO move to the case when Simulation is set to constant and use as default
|
||||
this->concentrations = MatrixXd::Constant(1, col, 20);
|
||||
this->alphaX = MatrixXd::Constant(1, col, 1);
|
||||
}
|
||||
@ -31,7 +30,6 @@ Grid::Grid(int row, int col) {
|
||||
this->deltaCol = double(this->domainCol)/double(this->col); // -> 1
|
||||
this->dim = 2;
|
||||
|
||||
// TODO move to the case when Simulation is set to constant and use as default
|
||||
this->concentrations = MatrixXd::Constant(row, col, 20);
|
||||
this->alphaX = MatrixXd::Constant(row, col, 1);
|
||||
this->alphaY = MatrixXd::Constant(row, col, 1);
|
||||
@ -145,6 +143,14 @@ void Grid::setDomain(int domainRow, int domainCol) {
|
||||
this->deltaCol = double(this->domainCol)/double(this->col);
|
||||
}
|
||||
|
||||
double Grid::getDelta() {
|
||||
if (dim != 1) {
|
||||
throw_invalid_argument("Grid is not one dimensional, you should probably use the 2D delta getters");
|
||||
}
|
||||
|
||||
return this->deltaCol;
|
||||
}
|
||||
|
||||
double Grid::getDeltaCol() {
|
||||
return this->deltaCol;
|
||||
}
|
||||
|
||||
45
src/README.md
Normal file
45
src/README.md
Normal file
@ -0,0 +1,45 @@
|
||||
<h1>src-Directory</h1>
|
||||
This is the src-directory that holds the source code to the implementation of the TUG framework.
|
||||
|
||||
<pre>
|
||||
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>
|
||||
|
||||
**src/** Directory with the source code.
|
||||
|
||||
**CMakeFiles/** Automatically generated directory by CMake.
|
||||
|
||||
**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.
|
||||
|
||||
**FTCS.cpp** Implementation of FTCS solution to heterogeneous diffusion in 1D and 2D.
|
||||
|
||||
**Grid.cpp** Implementation of Grid class, that holds all of the concentrations alpha coefficient in x- and y-direction.
|
||||
|
||||
**README.md** <i>This</i> file.
|
||||
|
||||
**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.
|
||||
|
||||
@ -19,14 +19,6 @@ Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid)
|
||||
this->timestep = -1; // error per default
|
||||
this->iterations = -1;
|
||||
this->innerIterations = 1;
|
||||
|
||||
// MDL no: we need to distinguish between "required dt" and
|
||||
// "number of (outer) iterations" at which the user needs an
|
||||
// output and the actual CFL-allowed timestep and consequently the
|
||||
// number of "inner" iterations which the explicit FTCS needs to
|
||||
// reach them. The following, at least at the moment, cannot be
|
||||
// computed here since "timestep" is not yet set when this
|
||||
// function is called. I brought everything into "FTCS_2D"!
|
||||
|
||||
this->csv_output = CSV_OUTPUT_OFF;
|
||||
this->console_output = CONSOLE_OUTPUT_OFF;
|
||||
@ -35,7 +27,6 @@ Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid)
|
||||
|
||||
void Simulation::setOutputCSV(CSV_OUTPUT csv_output) {
|
||||
if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) {
|
||||
// throw invalid_argument("Invalid CSV output option given!");
|
||||
throw_invalid_argument("Invalid CSV output option given!");
|
||||
}
|
||||
|
||||
@ -90,7 +81,7 @@ void Simulation::setTimestep(double timestep) {
|
||||
double CFL_Wiki = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia
|
||||
|
||||
cout << "FTCS_2D :: CFL condition MDL: " << CFL_MDL << endl;
|
||||
cout << "FTCS_2D :: CFL condition Wiki: " << CFL_Wiki << endl;
|
||||
// cout << "FTCS_2D :: CFL condition Wiki: " << CFL_Wiki << endl;
|
||||
cout << "FTCS_2D :: required dt=" << timestep << endl;
|
||||
|
||||
if (timestep > CFL_MDL) {
|
||||
@ -139,7 +130,6 @@ string Simulation::createCSVfile() {
|
||||
int appendIdent = 0;
|
||||
string appendIdentString;
|
||||
|
||||
// APPROACH_ROW_COL_ITERATIONS
|
||||
string approachString = (approach == 0) ? "FTCS" : "BTCS";
|
||||
string row = to_string(grid.getRow());
|
||||
string col = to_string(grid.getCol());
|
||||
@ -150,7 +140,7 @@ string Simulation::createCSVfile() {
|
||||
while (filesystem::exists(filename)) {
|
||||
appendIdent += 1;
|
||||
appendIdentString = to_string(appendIdent);
|
||||
filename = filename = approachString + "_" + row + "_" + col + "_" + numIterations + "-" + appendIdentString + ".csv";
|
||||
filename = approachString + "_" + row + "_" + col + "_" + numIterations + "-" + appendIdentString + ".csv";
|
||||
}
|
||||
|
||||
file.open(filename);
|
||||
@ -158,19 +148,14 @@ string Simulation::createCSVfile() {
|
||||
exit(1);
|
||||
}
|
||||
|
||||
// adds lines at the beginning of verbose output csv that represent the boundary conditions and their values
|
||||
// -1 in case of closed boundary
|
||||
if (csv_output == CSV_OUTPUT_XTREME) {
|
||||
//rows
|
||||
//cols
|
||||
//iterations
|
||||
//boundary left
|
||||
//boundary right
|
||||
//boundary top
|
||||
//boundary bottom
|
||||
file << row << endl;
|
||||
file << col << endl;
|
||||
file << numIterations << endl;
|
||||
// TODO
|
||||
// file << to_string(bc.printBoundarySide) << endl;
|
||||
IOFormat one_row(StreamPrecision, DontAlignCols, "", " ");
|
||||
file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) << endl; // boundary left
|
||||
file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) << endl; // boundary right
|
||||
file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) << endl; // boundary top
|
||||
file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) << endl; // boundary bottom
|
||||
file << endl << endl;
|
||||
}
|
||||
|
||||
@ -209,27 +194,26 @@ void Simulation::run() {
|
||||
filename = createCSVfile();
|
||||
}
|
||||
|
||||
auto begin = std::chrono::high_resolution_clock::now();
|
||||
|
||||
if (approach == FTCS_APPROACH) {
|
||||
auto begin = std::chrono::high_resolution_clock::now();
|
||||
progressbar bar(iterations * innerIterations);
|
||||
for (int i = 0; i < iterations * innerIterations; i++) {
|
||||
// MDL: distinguish between "outer" and "inner" iterations
|
||||
// std::cout << ":: run(): Outer iteration " << i+1 << "/" << iterations << endl;
|
||||
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
|
||||
printConcentrationsConsole();
|
||||
}
|
||||
if (csv_output == CSV_OUTPUT_VERBOSE) {
|
||||
if (csv_output >= CSV_OUTPUT_VERBOSE) {
|
||||
printConcentrationsCSV(filename);
|
||||
}
|
||||
|
||||
FTCS(this->grid, this->bc, this->timestep);
|
||||
bar.update();
|
||||
|
||||
if (i % (iterations * innerIterations / 100) == 0) {
|
||||
double percentage = (double)i / ((double)iterations * (double)innerIterations) * 100;
|
||||
if ((int)percentage % 10 == 0) {
|
||||
cout << "Progress: " << percentage << "%" << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
auto end = std::chrono::high_resolution_clock::now();
|
||||
auto milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
|
||||
|
||||
// MDL: meaningful stdout messages
|
||||
std::cout << ":: run() finished in " << milliseconds.count() << "ms" << endl;
|
||||
|
||||
} else if (approach == BTCS_APPROACH) {
|
||||
|
||||
@ -237,7 +221,7 @@ void Simulation::run() {
|
||||
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
|
||||
printConcentrationsConsole();
|
||||
}
|
||||
if (csv_output == CSV_OUTPUT_VERBOSE && i > 0) {
|
||||
if (csv_output >= CSV_OUTPUT_VERBOSE) {
|
||||
printConcentrationsCSV(filename);
|
||||
}
|
||||
|
||||
@ -247,12 +231,19 @@ void Simulation::run() {
|
||||
|
||||
}
|
||||
|
||||
auto end = std::chrono::high_resolution_clock::now();
|
||||
auto milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
|
||||
|
||||
if (this->console_output > CONSOLE_OUTPUT_OFF) {
|
||||
printConcentrationsConsole();
|
||||
}
|
||||
if (this->csv_output > CSV_OUTPUT_OFF) {
|
||||
printConcentrationsCSV(filename);
|
||||
}
|
||||
|
||||
if (this->time_measure > TIME_MEASURE_OFF) {
|
||||
string approachString = (approach == 0) ? "FTCS" : "BTCS";
|
||||
string dimString = (grid.getDim() == 1) ? "-1D" : "-2D";
|
||||
cout << approachString << dimString << ":: run() finished in " << milliseconds.count() << "ms" << endl;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -43,5 +43,5 @@ bool checkSimilarityV2(MatrixXd a, MatrixXd b, double maxDiff) {
|
||||
|
||||
MatrixXd diff = a - b;
|
||||
double maxCoeff = diff.maxCoeff();
|
||||
return maxCoeff < maxDiff;
|
||||
return abs(maxCoeff) < maxDiff;
|
||||
}
|
||||
@ -7,7 +7,7 @@
|
||||
// include the configured header file
|
||||
#include <testSimulation.hpp>
|
||||
|
||||
static Grid setupSimulation() {
|
||||
static Grid setupSimulation(APPROACH approach, double timestep, int iterations) {
|
||||
int row = 11;
|
||||
int col = 11;
|
||||
int domain_row = 10;
|
||||
@ -46,9 +46,10 @@ static Grid setupSimulation() {
|
||||
|
||||
|
||||
// Simulation
|
||||
Simulation sim = Simulation(grid, bc, FTCS_APPROACH);
|
||||
sim.setTimestep(0.001);
|
||||
sim.setIterations(7000);
|
||||
Simulation sim = Simulation(grid, bc, approach);
|
||||
sim.setOutputConsole(CONSOLE_OUTPUT_ON);
|
||||
sim.setTimestep(timestep);
|
||||
sim.setIterations(iterations);
|
||||
sim.run();
|
||||
|
||||
// RUN
|
||||
@ -56,21 +57,29 @@ static Grid setupSimulation() {
|
||||
|
||||
}
|
||||
|
||||
TEST_CASE("equality to reference matrix") {
|
||||
TEST_CASE("equality to reference matrix with FTCS") {
|
||||
// set string from the header file
|
||||
string test_path = testSimulationCSVDir;
|
||||
MatrixXd reference = CSV2Eigen(test_path);
|
||||
Grid grid = setupSimulation();
|
||||
Grid grid = setupSimulation(FTCS_APPROACH, 0.001, 7000);
|
||||
CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true);
|
||||
}
|
||||
|
||||
TEST_CASE("equality to reference matrix with BTCS") {
|
||||
// set string from the header file
|
||||
string test_path = testSimulationCSVDir;
|
||||
MatrixXd reference = CSV2Eigen(test_path);
|
||||
Grid grid = setupSimulation(BTCS_APPROACH, 1, 7);
|
||||
CHECK(checkSimilarityV2(reference, grid.getConcentrations(), 0.01) == true);
|
||||
}
|
||||
|
||||
TEST_CASE("Initialize environment"){
|
||||
int rc = 12;
|
||||
Grid grid(rc, rc);
|
||||
Boundary boundary(grid);
|
||||
|
||||
CHECK_NOTHROW(Simulation sim(grid, boundary, FTCS_APPROACH));
|
||||
}
|
||||
CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH));
|
||||
}
|
||||
|
||||
TEST_CASE("Simulation environment"){
|
||||
int rc = 12;
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user