mirror of
https://git.gfz-potsdam.de/naaice/tug.git
synced 2025-12-13 09:28:23 +01:00
Merge branch 'dynamic-boundary-conditions-functionality' into 'hannes-philipp'
Dynamic boundary conditions functionality See merge request naaice/tug!6
This commit is contained in:
Hannes Martin Signer
commit
2d293469ff
@ -7,7 +7,7 @@ set(CMAKE_CXX_STANDARD 17)
|
||||
|
||||
find_package(Eigen3 REQUIRED NO_MODULE)
|
||||
find_package(OpenMP)
|
||||
find_package(easy_profiler REQUIRED)
|
||||
find_package(easy_profiler)
|
||||
|
||||
## SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -mfma")
|
||||
option(TUG_USE_OPENMP "Compile with OpenMP support" ON)
|
||||
|
||||
@ -9,4 +9,4 @@ target_link_libraries(second_example tug)
|
||||
target_link_libraries(boundary_example1D tug)
|
||||
target_link_libraries(FTCS_2D_proto_example tug)
|
||||
target_link_libraries(FTCS_1D_proto_example tug)
|
||||
target_link_libraries(FTCS_2D_proto_example easy_profiler)
|
||||
# target_link_libraries(FTCS_2D_proto_example easy_profiler)
|
||||
@ -1,3 +1,4 @@
|
||||
#include "tug/Boundary.hpp"
|
||||
#include <tug/Simulation.hpp>
|
||||
|
||||
int main(int argc, char *argv[]) {
|
||||
@ -18,7 +19,9 @@ int main(int argc, char *argv[]) {
|
||||
// ******************
|
||||
|
||||
// create a boundary with constant values
|
||||
Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
|
||||
Boundary bc = Boundary(grid);
|
||||
bc.setBoundarySideConstant(BC_SIDE_LEFT, 1);
|
||||
bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1);
|
||||
|
||||
|
||||
// ************************
|
||||
|
||||
@ -7,17 +7,17 @@
|
||||
*/
|
||||
|
||||
#include <tug/Simulation.hpp>
|
||||
#include <easy/profiler.h>
|
||||
#define EASY_PROFILER_ENABLE ::profiler::setEnabled(true);
|
||||
// #include <easy/profiler.h>
|
||||
// #define EASY_PROFILER_ENABLE ::profiler::setEnabled(true);
|
||||
|
||||
|
||||
int main(int argc, char *argv[]) {
|
||||
EASY_PROFILER_ENABLE;
|
||||
profiler::startListen();
|
||||
// EASY_PROFILER_ENABLE;
|
||||
// profiler::startListen();
|
||||
// **************
|
||||
// **** GRID ****
|
||||
// **************
|
||||
profiler::startListen();
|
||||
// profiler::startListen();
|
||||
// create a grid with a 20 x 20 field
|
||||
int row = 20;
|
||||
int col = 20;
|
||||
@ -29,8 +29,8 @@ int main(int argc, char *argv[]) {
|
||||
// (optional) set the concentrations, e.g.:
|
||||
// MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value
|
||||
// grid.setConcentrations(concentrations);
|
||||
MatrixXd concentrations = MatrixXd::Constant(row, col,1);
|
||||
concentrations(0,0) = 2000;
|
||||
MatrixXd concentrations = MatrixXd::Constant(row,col,0);
|
||||
concentrations(0,0) = 1999;
|
||||
grid.setConcentrations(concentrations);
|
||||
|
||||
|
||||
@ -45,7 +45,12 @@ int main(int argc, char *argv[]) {
|
||||
// ******************
|
||||
|
||||
// create a boundary with constant values
|
||||
Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
|
||||
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);
|
||||
|
||||
|
||||
// (optional) set boundary condition values for one side, e.g.:
|
||||
// VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value
|
||||
@ -69,7 +74,7 @@ int main(int argc, char *argv[]) {
|
||||
simulation.setTimestep(0.1); // timestep
|
||||
|
||||
// (optional) set the number of iterations
|
||||
simulation.setIterations(2);
|
||||
simulation.setIterations(10000);
|
||||
|
||||
// (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
|
||||
simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
|
||||
@ -78,9 +83,9 @@ int main(int argc, char *argv[]) {
|
||||
|
||||
// run the simulation
|
||||
|
||||
EASY_BLOCK("SIMULATION")
|
||||
// EASY_BLOCK("SIMULATION")
|
||||
simulation.run();
|
||||
EASY_END_BLOCK;
|
||||
profiler::dumpBlocksToFile("test_profile.prof");
|
||||
profiler::stopListen();
|
||||
// EASY_END_BLOCK;
|
||||
// profiler::dumpBlocksToFile("test_profile.prof");
|
||||
// profiler::stopListen();
|
||||
}
|
||||
@ -2,6 +2,7 @@
|
||||
#define BOUNDARY_H_
|
||||
|
||||
#include <Eigen/Core>
|
||||
#include <cstddef>
|
||||
#include "Grid.hpp"
|
||||
|
||||
using namespace std;
|
||||
@ -19,26 +20,21 @@ enum BC_SIDE {
|
||||
BC_SIDE_BOTTOM
|
||||
};
|
||||
|
||||
/*******************************************
|
||||
class WallElement {
|
||||
class BoundaryElement {
|
||||
public:
|
||||
WallElement() {
|
||||
this->type = BC_TYPE_CLOSED;
|
||||
this->value = 0;
|
||||
}
|
||||
// bc type closed
|
||||
BoundaryElement();
|
||||
|
||||
WallElement(double value) {
|
||||
this->type = BC_TYPE_CONSTANT;
|
||||
this->value = value;
|
||||
}
|
||||
// bc type constant
|
||||
BoundaryElement(double value);
|
||||
|
||||
BC_TYPE getType() {
|
||||
return this->type;
|
||||
}
|
||||
void setType(BC_TYPE type);
|
||||
|
||||
double getValue() {
|
||||
return this->value;
|
||||
}
|
||||
void setValue(double value);
|
||||
|
||||
BC_TYPE getType();
|
||||
|
||||
double getValue();
|
||||
|
||||
private:
|
||||
BC_TYPE type;
|
||||
@ -47,16 +43,22 @@ class WallElement {
|
||||
|
||||
class BoundaryWall {
|
||||
public:
|
||||
BoundaryWall(int length) {
|
||||
// create array with length many wall elements
|
||||
}
|
||||
BoundaryWall(int length);
|
||||
|
||||
void setWall(BC_TYPE type, double value = NAN);
|
||||
|
||||
vector<BoundaryElement> getWall();
|
||||
|
||||
void setBoundaryElement(int index, BC_TYPE type, double value = NAN);
|
||||
|
||||
BoundaryElement getBoundaryElement();
|
||||
|
||||
private:
|
||||
BC_SIDE side;
|
||||
int length;
|
||||
vector<WallElement> wall;
|
||||
vector<BoundaryElement> wall;
|
||||
|
||||
};
|
||||
***********************/
|
||||
|
||||
class Boundary {
|
||||
public:
|
||||
@ -67,40 +69,29 @@ class Boundary {
|
||||
* @param grid
|
||||
* @param type
|
||||
*/
|
||||
Boundary(Grid grid, BC_TYPE type);
|
||||
Boundary(Grid grid);
|
||||
|
||||
/**
|
||||
* @brief Get the Boundary Condition Type object
|
||||
*
|
||||
* @return auto
|
||||
*/
|
||||
BC_TYPE getBoundaryConditionType();
|
||||
void setBoundarySideClosed(BC_SIDE side);
|
||||
|
||||
/**
|
||||
* @brief Set the Boundary Condition Value object
|
||||
*
|
||||
* @param side
|
||||
* @param values
|
||||
*/
|
||||
void setBoundaryConditionValue(BC_SIDE side, VectorXd values);
|
||||
void setBoundarySideConstant(BC_SIDE side, double value);
|
||||
|
||||
/**
|
||||
* @brief Get the Boundary Condition Value object
|
||||
*
|
||||
* @param side
|
||||
* @return auto
|
||||
*/
|
||||
VectorXd getBoundaryConditionValue(BC_SIDE side);
|
||||
void setBoundaryElementClosed(BC_SIDE side, int index);
|
||||
|
||||
void setBoundaryElementConstant(BC_SIDE side, int index, double value);
|
||||
|
||||
vector<BoundaryElement> getBoundarySide(BC_SIDE side);
|
||||
|
||||
BoundaryElement getBoundaryElement(BC_SIDE side, int index);
|
||||
|
||||
BC_TYPE getBoundaryElementType(BC_SIDE side, int index);
|
||||
|
||||
double getBoundaryElementValue(BC_SIDE side, int index);
|
||||
|
||||
private:
|
||||
Grid grid;
|
||||
|
||||
// need a way to save the bc type and value for each single 'boundary cell'
|
||||
// perhaps an array for each side with structs containing the bc type as well as a value
|
||||
// or another object that contains one boundary side
|
||||
BC_TYPE type;
|
||||
VectorXd left, right, top, bottom;
|
||||
|
||||
vector<vector<BoundaryElement>> boundaries;
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
123
src/Boundary.cpp
123
src/Boundary.cpp
@ -1,68 +1,85 @@
|
||||
#include "tug/BoundaryCondition.hpp"
|
||||
#include <iostream>
|
||||
#include <omp.h>
|
||||
#include <tug/Boundary.hpp>
|
||||
#include <stdexcept>
|
||||
|
||||
using namespace std;
|
||||
|
||||
Boundary::Boundary(Grid grid, BC_TYPE type) : grid(grid) {
|
||||
//probably to DEBUG assignment grid
|
||||
BoundaryElement::BoundaryElement() {
|
||||
this->type = BC_TYPE_CLOSED;
|
||||
this->value = NAN;
|
||||
}
|
||||
|
||||
BoundaryElement::BoundaryElement(double value) {
|
||||
this->type = BC_TYPE_CONSTANT;
|
||||
this->value = value;
|
||||
}
|
||||
|
||||
void BoundaryElement::setType(BC_TYPE type) {
|
||||
this->type = type;
|
||||
|
||||
if (type == BC_TYPE_CONSTANT) {
|
||||
if (grid.getDim() == 1) {
|
||||
this->left = VectorXd::Constant(1, 1);
|
||||
this->right = VectorXd::Constant(1, 1);
|
||||
} else if (grid.getDim() == 2) {
|
||||
this->left = VectorXd::Constant(grid.getRow(), 1);
|
||||
this->right = VectorXd::Constant(grid.getRow(), 1);
|
||||
this->top = VectorXd::Constant(grid.getCol(), 1);
|
||||
this->bottom = VectorXd::Constant(grid.getCol(), 1);
|
||||
|
||||
} else {
|
||||
throw invalid_argument("Dimension must be 1 or 2!");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
BC_TYPE Boundary::getBoundaryConditionType() {
|
||||
return this->type;
|
||||
void BoundaryElement::setValue(double value) {
|
||||
this->value = value;
|
||||
}
|
||||
|
||||
void Boundary::setBoundaryConditionValue(BC_SIDE side, VectorXd values) {
|
||||
if (type != BC_TYPE_CONSTANT) {
|
||||
// TODO check if correct way for handling warning
|
||||
cerr << "Values will not be used, wrong BC_TYPE!";
|
||||
}
|
||||
BC_TYPE BoundaryElement::getType() {
|
||||
return this->type;
|
||||
}
|
||||
|
||||
switch (side) {
|
||||
case BC_SIDE_LEFT:
|
||||
this->left = values;
|
||||
break;
|
||||
case BC_SIDE_RIGHT:
|
||||
this->right = values;
|
||||
break;
|
||||
case BC_SIDE_TOP:
|
||||
this->top = values;
|
||||
break;
|
||||
case BC_SIDE_BOTTOM:
|
||||
this->bottom = values;
|
||||
break;
|
||||
default:
|
||||
throw invalid_argument("Invalid side given!");
|
||||
double BoundaryElement::getValue() {
|
||||
return this->value;
|
||||
}
|
||||
|
||||
Boundary::Boundary(Grid grid) : grid(grid) {
|
||||
//probably to DEBUG assignment grid
|
||||
|
||||
if (grid.getDim() == 1) {
|
||||
this->boundaries = vector<vector<BoundaryElement>>(2);
|
||||
|
||||
this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement());
|
||||
this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement());
|
||||
} else if (grid.getDim() == 2) {
|
||||
this->boundaries = vector<vector<BoundaryElement>>(4);
|
||||
|
||||
this->boundaries[BC_SIDE_LEFT] = vector<BoundaryElement>(grid.getRow(), BoundaryElement());
|
||||
this->boundaries[BC_SIDE_RIGHT] = vector<BoundaryElement>(grid.getRow(), BoundaryElement());
|
||||
this->boundaries[BC_SIDE_TOP] = vector<BoundaryElement>(grid.getCol(), BoundaryElement());
|
||||
this->boundaries[BC_SIDE_BOTTOM] = vector<BoundaryElement>(grid.getCol(), BoundaryElement());
|
||||
}
|
||||
}
|
||||
|
||||
VectorXd Boundary::getBoundaryConditionValue(BC_SIDE side) {
|
||||
switch (side) {
|
||||
case BC_SIDE_LEFT:
|
||||
return this->left;
|
||||
case BC_SIDE_RIGHT:
|
||||
return this->right;
|
||||
case BC_SIDE_TOP:
|
||||
return this->top;
|
||||
case BC_SIDE_BOTTOM:
|
||||
return this->bottom;
|
||||
default:
|
||||
throw invalid_argument("Invalid side given!");
|
||||
}
|
||||
}
|
||||
void Boundary::setBoundarySideClosed(BC_SIDE side) {
|
||||
this->boundaries[side] = vector<BoundaryElement>(grid.getRow(), BoundaryElement());
|
||||
}
|
||||
|
||||
void Boundary::setBoundarySideConstant(BC_SIDE side, double value) {
|
||||
this->boundaries[side] = vector<BoundaryElement>(grid.getRow(), BoundaryElement(value));
|
||||
}
|
||||
|
||||
void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) {
|
||||
this->boundaries[side][index].setType(BC_TYPE_CLOSED);
|
||||
}
|
||||
|
||||
void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, double value) {
|
||||
this->boundaries[side][index].setType(BC_TYPE_CONSTANT);
|
||||
this->boundaries[side][index].setValue(value);
|
||||
}
|
||||
|
||||
vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
|
||||
return this->boundaries[side];
|
||||
}
|
||||
|
||||
BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) {
|
||||
return this->boundaries[side][index];
|
||||
}
|
||||
|
||||
BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) {
|
||||
return this->boundaries[side][index].getType();
|
||||
}
|
||||
|
||||
double Boundary::getBoundaryElementValue(BC_SIDE side, int index) {
|
||||
return this->boundaries[side][index].getValue();
|
||||
}
|
||||
|
||||
|
||||
138
src/FTCS.cpp
138
src/FTCS.cpp
@ -1,9 +1,11 @@
|
||||
#include "TugUtils.hpp"
|
||||
#include <cstddef>
|
||||
#include <tug/Boundary.hpp>
|
||||
#include <iostream>
|
||||
|
||||
using namespace std;
|
||||
|
||||
|
||||
double calcAlphaIntercell(double alpha1, double alpha2, bool useHarmonic = false) {
|
||||
if (useHarmonic) {
|
||||
return 2 / ((1/alpha1) + (1/alpha2));
|
||||
@ -57,21 +59,37 @@ double calcHorizontalChangeLeftBoundaryConstant(Grid grid, Boundary bc, int row,
|
||||
+ 2 * grid.getAlphaX()(row,col)
|
||||
)
|
||||
* grid.getConcentrations()(row,col)
|
||||
+ 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row);
|
||||
+ 2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_LEFT, row);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
double calcHorizontalChangeLeftBoundaryClosed() {
|
||||
return 0;
|
||||
double calcHorizontalChangeLeftBoundaryClosed(Grid grid, int row, int col) {
|
||||
|
||||
double result =
|
||||
calcAlphaIntercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col))
|
||||
* (grid.getConcentrations()(row, col+1) - grid.getConcentrations()(row, col));
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
double calcHorizontalChangeLeftBoundary(Grid grid, Boundary bc, int row, int col) {
|
||||
if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) {
|
||||
return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col);
|
||||
} else if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CLOSED) {
|
||||
return calcHorizontalChangeLeftBoundaryClosed(grid, row, col);
|
||||
} else {
|
||||
throw_invalid_argument("Undefined Boundary Condition Type!");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double calcHorizontalChangeRightBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
|
||||
|
||||
double result =
|
||||
2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
|
||||
2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT, row)
|
||||
- (
|
||||
calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
|
||||
+ 2 * grid.getAlphaX()(row,col)
|
||||
@ -84,8 +102,24 @@ double calcHorizontalChangeRightBoundaryConstant(Grid grid, Boundary bc, int row
|
||||
}
|
||||
|
||||
|
||||
double calcHorizontalChangeRightBoundaryClosed() {
|
||||
return 0;
|
||||
double calcHorizontalChangeRightBoundaryClosed(Grid grid, int row, int col) {
|
||||
|
||||
double result =
|
||||
- (calcAlphaIntercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col))
|
||||
* (grid.getConcentrations()(row, col) - grid.getConcentrations()(row, col-1)));
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
double calcHorizontalChangeRightBoundary(Grid grid, Boundary bc, int row, int col) {
|
||||
if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) {
|
||||
return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col);
|
||||
} else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CLOSED) {
|
||||
return calcHorizontalChangeRightBoundaryClosed(grid, row, col);
|
||||
} else {
|
||||
throw_invalid_argument("Undefined Boundary Condition Type!");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -99,21 +133,37 @@ double calcVerticalChangeTopBoundaryConstant(Grid grid, Boundary bc, int row, in
|
||||
+ 2 * grid.getAlphaY()(row, col)
|
||||
)
|
||||
* grid.getConcentrations()(row, col)
|
||||
+ 2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col);
|
||||
+ 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_TOP, col);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
double calcVerticalChangeTopBoundaryClosed() {
|
||||
return 0;
|
||||
double calcVerticalChangeTopBoundaryClosed(Grid grid, int row, int col) {
|
||||
|
||||
double result =
|
||||
calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getConcentrations()(row, col))
|
||||
* (grid.getConcentrations()(row+1, col) - grid.getConcentrations()(row, col));
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
double calcVerticalChangeTopBoundary(Grid grid, Boundary bc, int row, int col) {
|
||||
if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) {
|
||||
return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col);
|
||||
} else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) {
|
||||
return calcVerticalChangeTopBoundaryClosed(grid, row, col);
|
||||
} else {
|
||||
throw_invalid_argument("Undefined Boundary Condition Type!");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double calcVerticalChangeBottomBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
|
||||
|
||||
double result =
|
||||
2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
|
||||
2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col)
|
||||
- (
|
||||
calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col))
|
||||
+ 2 * grid.getAlphaY()(row, col)
|
||||
@ -126,8 +176,24 @@ double calcVerticalChangeBottomBoundaryConstant(Grid grid, Boundary bc, int row,
|
||||
}
|
||||
|
||||
|
||||
double calcVerticalChangeBottomBoundaryClosed() {
|
||||
return 0;
|
||||
double calcVerticalChangeBottomBoundaryClosed(Grid grid, int row, int col) {
|
||||
|
||||
double result =
|
||||
- (calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col))
|
||||
* (grid.getConcentrations()(row, col) - grid.getConcentrations()(row-1, col)));
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
double calcVerticalChangeBottomBoundary(Grid grid, Boundary bc, int row, int col) {
|
||||
if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) {
|
||||
return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col);
|
||||
} else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) {
|
||||
return calcVerticalChangeBottomBoundaryClosed(grid, row, col);
|
||||
} else {
|
||||
throw_invalid_argument("Undefined Boundary Condition Type!");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -155,7 +221,7 @@ MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) {
|
||||
concentrations_t1(row, col) = grid.getConcentrations()(row,col)
|
||||
+ timestep / (deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
|
||||
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
@ -165,7 +231,7 @@ MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) {
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ timestep / (deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
|
||||
calcHorizontalChangeRightBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
@ -207,7 +273,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
|
||||
concentrations_t1(row, col) = grid.getConcentrations()(row,col)
|
||||
+ timestep / (deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
|
||||
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep / (deltaRow*deltaRow)
|
||||
* (
|
||||
@ -222,7 +288,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ timestep / (deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
|
||||
calcHorizontalChangeRightBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep / (deltaRow*deltaRow)
|
||||
* (
|
||||
@ -238,7 +304,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
|
||||
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
|
||||
+ timestep / (deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
|
||||
calcVerticalChangeTopBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep / (deltaCol*deltaCol)
|
||||
* (
|
||||
@ -253,7 +319,7 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
|
||||
concentrations_t1(row, col) = grid.getConcentrations()(row, col)
|
||||
+ timestep / (deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
|
||||
calcVerticalChangeBottomBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep / (deltaCol*deltaCol)
|
||||
* (
|
||||
@ -268,11 +334,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ timestep/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
|
||||
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
|
||||
calcVerticalChangeTopBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
@ -282,11 +348,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ timestep/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
|
||||
calcHorizontalChangeRightBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
|
||||
calcVerticalChangeTopBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
@ -296,11 +362,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ timestep/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
|
||||
calcHorizontalChangeLeftBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
|
||||
calcVerticalChangeBottomBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
@ -310,11 +376,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
|
||||
concentrations_t1(row,col) = grid.getConcentrations()(row,col)
|
||||
+ timestep/(deltaCol*deltaCol)
|
||||
* (
|
||||
calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
|
||||
calcHorizontalChangeRightBoundary(grid, bc, row, col)
|
||||
)
|
||||
+ timestep/(deltaRow*deltaRow)
|
||||
* (
|
||||
calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
|
||||
calcVerticalChangeBottomBoundary(grid, bc, row, col)
|
||||
)
|
||||
;
|
||||
|
||||
@ -324,29 +390,11 @@ MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
|
||||
|
||||
|
||||
MatrixXd FTCS(Grid grid, Boundary bc, double timestep) {
|
||||
// inner cells
|
||||
// TODO only the boundary cells are different in constant and closed case
|
||||
|
||||
// if 1D:
|
||||
// do inner cells
|
||||
// do left boundary according to bc type
|
||||
// do right boundary according to bc type
|
||||
// if 2D:
|
||||
// do inner cells
|
||||
// do left boundaries according to bc type
|
||||
// do right boundaries according to bc type
|
||||
// ...
|
||||
|
||||
|
||||
if (grid.getDim() == 1) {
|
||||
return FTCS_1D(grid, bc, timestep);
|
||||
} else {
|
||||
return FTCS_2D(grid, bc, timestep);
|
||||
}
|
||||
|
||||
// checking the boundary condition type first does not work
|
||||
// if the boundary condition types change dynamically for a grid
|
||||
// meaning:
|
||||
// - boundary condition type needs to be checked for every single boundary cell
|
||||
// -> this check is last in order
|
||||
// -
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user