diff --git a/examples/BTCS_1D_proto_example.cpp b/examples/BTCS_1D_proto_example.cpp index fb04966..f7e7a18 100644 --- a/examples/BTCS_1D_proto_example.cpp +++ b/examples/BTCS_1D_proto_example.cpp @@ -1,48 +1,48 @@ #include int main(int argc, char *argv[]) { - // ************** - // **** GRID **** - // ************** + // ************** + // **** GRID **** + // ************** - // create a linear grid with 20 cells - int cells = 20; - Grid grid = Grid(cells); + // create a linear grid with 20 cells + int cells = 20; + Grid grid = Grid(cells); - MatrixXd concentrations = MatrixXd::Constant(1,20,0); - concentrations(0,0) = 2000; - // TODO add option to set concentrations with a vector in 1D case - grid.setConcentrations(concentrations); + MatrixXd concentrations = MatrixXd::Constant(1, 20, 0); + concentrations(0, 0) = 2000; + // TODO add option to set concentrations with a vector in 1D case + grid.setConcentrations(concentrations); + // ****************** + // **** BOUNDARY **** + // ****************** - // ****************** - // **** BOUNDARY **** - // ****************** + // create a boundary with constant values + Boundary bc = Boundary(grid); + bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); + bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); - // create a boundary with constant values - Boundary bc = Boundary(grid); - bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); - bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + // ************************ + // **** SIMULATION ENV **** + // ************************ + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach - // ************************ - // **** SIMULATION ENV **** - // ************************ + // set the timestep of the simulation + simulation.setTimestep(0.1); // timestep - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + // set the number of iterations + simulation.setIterations(100); - // set the timestep of the simulation - simulation.setTimestep(0.1); // timestep + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - // set the number of iterations - simulation.setIterations(100); + // **** RUN SIMULATION **** - // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - - // **** RUN SIMULATION **** - - // run the simulation - simulation.run(); + // run the simulation + simulation.run(); } \ No newline at end of file diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index c7b1a17..f65f572 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -1,84 +1,82 @@ #include int main(int argc, char *argv[]) { - // EASY_PROFILER_ENABLE; - // profiler::startListen(); - // ************** - // **** GRID **** - // ************** - // profiler::startListen(); - // create a grid with a 20 x 20 field - int row = 40; - int col = 50; - Grid grid = Grid(row,col); + // EASY_PROFILER_ENABLE; + // profiler::startListen(); + // ************** + // **** GRID **** + // ************** + // profiler::startListen(); + // create a grid with a 20 x 20 field + int row = 40; + int col = 50; + Grid grid = Grid(row, col); - // (optional) set the domain, e.g.: - // grid.setDomain(20, 20); + // (optional) set the domain, e.g.: + // grid.setDomain(20, 20); - // (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,0); - concentrations(10,10) = 2000; - grid.setConcentrations(concentrations); - + // (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, 0); + concentrations(10, 10) = 2000; + grid.setConcentrations(concentrations); - // (optional) set alphas of the grid, e.g.: - // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value - // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value - // grid.setAlpha(alphax, alphay); + // (optional) set alphas of the grid, e.g.: + // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value + // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value + // grid.setAlpha(alphax, alphay); + // ****************** + // **** BOUNDARY **** + // ****************** - // ****************** - // **** BOUNDARY **** - // ****************** + // create a boundary with constant values + Boundary bc = Boundary(grid); + 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); - // create a boundary with constant values - Boundary bc = Boundary(grid); - 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.: + // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value + // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values + // VectorXd bc_zero_values = VectorXd::Constant(20,0); + // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); + // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); + // VectorXd bc_front_values = VectorXd::Constant(20,2000); + // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); + // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); + // ************************ + // **** SIMULATION ENV **** + // ************************ - // (optional) set boundary condition values for one side, e.g.: - // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value - // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values - // VectorXd bc_zero_values = VectorXd::Constant(20,0); - // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); - // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); - // VectorXd bc_front_values = VectorXd::Constant(20,2000); - // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); - // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + // set the timestep of the simulation + simulation.setTimestep(0.1); // timestep - // ************************ - // **** SIMULATION ENV **** - // ************************ + // set the number of iterations + simulation.setIterations(300); - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_XTREME); - // set the timestep of the simulation - simulation.setTimestep(0.1); // timestep + // **** RUN SIMULATION **** - // set the number of iterations - simulation.setIterations(300); + // run the simulation - // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_XTREME); - - // **** RUN SIMULATION **** - - // run the simulation - - // EASY_BLOCK("SIMULATION") - simulation.run(); - // EASY_END_BLOCK; - // profiler::dumpBlocksToFile("test_profile.prof"); - // profiler::stopListen(); + // EASY_BLOCK("SIMULATION") + simulation.run(); + // EASY_END_BLOCK; + // profiler::dumpBlocksToFile("test_profile.prof"); + // profiler::stopListen(); } \ No newline at end of file diff --git a/examples/CRNI_2D_proto_example.cpp b/examples/CRNI_2D_proto_example.cpp index c976b51..770fd42 100644 --- a/examples/CRNI_2D_proto_example.cpp +++ b/examples/CRNI_2D_proto_example.cpp @@ -1,24 +1,24 @@ #include int main(int argc, char *argv[]) { - int row = 20; - int col = 20; - Grid grid(row, col); + int row = 20; + int col = 20; + Grid grid(row, col); - MatrixXd concentrations = MatrixXd::Constant(row,col,0); - concentrations(10,10) = 2000; - grid.setConcentrations(concentrations); + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(10, 10) = 2000; + grid.setConcentrations(concentrations); - Boundary bc = Boundary(grid); - bc.setBoundarySideClosed(BC_SIDE_LEFT); - bc.setBoundarySideClosed(BC_SIDE_RIGHT); - bc.setBoundarySideClosed(BC_SIDE_TOP); - bc.setBoundarySideClosed(BC_SIDE_BOTTOM); + Boundary bc = Boundary(grid); + bc.setBoundarySideClosed(BC_SIDE_LEFT); + bc.setBoundarySideClosed(BC_SIDE_RIGHT); + bc.setBoundarySideClosed(BC_SIDE_TOP); + bc.setBoundarySideClosed(BC_SIDE_BOTTOM); - Simulation simulation = Simulation(grid, bc, CRANK_NICOLSON_APPROACH); - simulation.setTimestep(0.1); - simulation.setIterations(50); - simulation.setOutputCSV(CSV_OUTPUT_XTREME); + Simulation simulation = Simulation(grid, bc, CRANK_NICOLSON_APPROACH); + simulation.setTimestep(0.1); + simulation.setIterations(50); + simulation.setOutputCSV(CSV_OUTPUT_XTREME); - simulation.run(); + simulation.run(); } \ No newline at end of file diff --git a/examples/FTCS_1D_proto_example.cpp b/examples/FTCS_1D_proto_example.cpp index 3ac0df6..0b2ee83 100644 --- a/examples/FTCS_1D_proto_example.cpp +++ b/examples/FTCS_1D_proto_example.cpp @@ -2,46 +2,46 @@ #include int main(int argc, char *argv[]) { - // ************** - // **** GRID **** - // ************** + // ************** + // **** GRID **** + // ************** - // create a linear grid with 20 cells - int cells = 20; - Grid grid = Grid(cells); + // create a linear grid with 20 cells + int cells = 20; + Grid grid = Grid(cells); - MatrixXd concentrations = MatrixXd::Constant(1,20,20); - grid.setConcentrations(concentrations); + MatrixXd concentrations = MatrixXd::Constant(1, 20, 20); + grid.setConcentrations(concentrations); + // ****************** + // **** BOUNDARY **** + // ****************** - // ****************** - // **** BOUNDARY **** - // ****************** + // create a boundary with constant values + Boundary bc = Boundary(grid); + bc.setBoundarySideConstant(BC_SIDE_LEFT, 1); + bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1); - // create a boundary with constant values - Boundary bc = Boundary(grid); - bc.setBoundarySideConstant(BC_SIDE_LEFT, 1); - bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1); + // ************************ + // **** SIMULATION ENV **** + // ************************ + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach - // ************************ - // **** SIMULATION ENV **** - // ************************ + // (optional) set the timestep of the simulation + simulation.setTimestep(0.1); // timestep - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + // (optional) set the number of iterations + simulation.setIterations(100); - // (optional) set the timestep of the simulation - simulation.setTimestep(0.1); // timestep + // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_OFF); - // (optional) set the number of iterations - simulation.setIterations(100); + // **** RUN SIMULATION **** - // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_OFF); - - // **** RUN SIMULATION **** - - // run the simulation - simulation.run(); + // run the simulation + simulation.run(); } \ No newline at end of file diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index 4ed9599..a297734 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -1,8 +1,9 @@ /** * @file FTCS_2D_proto_closed_mdl.cpp * @author Hannes Signer, Philipp Ungrund, MDL - * @brief Creates a TUG simulation in 2D with FTCS approach and closed boundary condition; optional command line argument: number of cols and rows - * + * @brief Creates a TUG simulation in 2D with FTCS approach and closed boundary + * condition; optional command line argument: number of cols and rows + * */ #include @@ -11,75 +12,76 @@ int main(int argc, char *argv[]) { - int row = 64; + int row = 64; - if (argc == 2) { - // no cmd line argument, take col=row=64 - row = atoi(argv[1]); - } - int col=row; + if (argc == 2) { + // no cmd line argument, take col=row=64 + row = atoi(argv[1]); + } + int col = row; - std::cout << "Nrow =" << row << std::endl; - // ************** - // **** GRID **** - // ************** + std::cout << "Nrow =" << row << std::endl; + // ************** + // **** GRID **** + // ************** - // create a grid with a 20 x 20 field - int n2 = row/2-1; - Grid grid = Grid(row,col); + // create a grid with a 20 x 20 field + int n2 = row / 2 - 1; + Grid grid = Grid(row, col); - // (optional) set the domain, e.g.: - // grid.setDomain(20, 20); + // (optional) set the domain, e.g.: + // grid.setDomain(20, 20); - // (optional) set the concentrations, e.g.: - // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value - MatrixXd concentrations = MatrixXd::Constant(row,col,0); - concentrations(n2,n2) = 1; - concentrations(n2,n2+1) = 1; - concentrations(n2+1,n2) = 1; - concentrations(n2+1,n2+1) = 1; - grid.setConcentrations(concentrations); + // (optional) set the concentrations, e.g.: + // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // + // #row,#col,value + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(n2, n2) = 1; + concentrations(n2, n2 + 1) = 1; + concentrations(n2 + 1, n2) = 1; + concentrations(n2 + 1, n2 + 1) = 1; + grid.setConcentrations(concentrations); - // (optional) set alphas of the grid, e.g.: - MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value - MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value - grid.setAlpha(alphax, alphay); + // (optional) set alphas of the grid, e.g.: + MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value + MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value + grid.setAlpha(alphax, alphay); + // ****************** + // **** BOUNDARY **** + // ****************** - // ****************** - // **** BOUNDARY **** - // ****************** + // create a boundary with constant values + Boundary bc = Boundary(grid); - // create a boundary with constant values - Boundary bc = Boundary(grid); + // (optional) set boundary condition values for one side, e.g.: + bc.setBoundarySideClosed(BC_SIDE_LEFT); // side,values + bc.setBoundarySideClosed(BC_SIDE_RIGHT); + bc.setBoundarySideClosed(BC_SIDE_TOP); + bc.setBoundarySideClosed(BC_SIDE_BOTTOM); - // (optional) set boundary condition values for one side, e.g.: - bc.setBoundarySideClosed(BC_SIDE_LEFT); // side,values - bc.setBoundarySideClosed(BC_SIDE_RIGHT); - bc.setBoundarySideClosed(BC_SIDE_TOP); - bc.setBoundarySideClosed(BC_SIDE_BOTTOM); + // ************************ + // **** SIMULATION ENV **** + // ************************ + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach - // ************************ - // **** SIMULATION ENV **** - // ************************ + // set the timestep of the simulation + simulation.setTimestep(10000); // timestep - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + // set the number of iterations + simulation.setIterations(100); - // set the timestep of the simulation - simulation.setTimestep(10000); // timestep + // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - // set the number of iterations - simulation.setIterations(100); + // **** RUN SIMULATION **** - // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - - // **** RUN SIMULATION **** + // run the simulation + simulation.run(); - // run the simulation - simulation.run(); - - return 0; + return 0; } diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index e65d425..058810d 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -1,92 +1,88 @@ /** * @file FTCS_2D_proto_example.cpp * @author Hannes Signer, Philipp Ungrund - * @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach - * and constant boundary condition - * + * @brief Creates a prototypical standard TUG simulation in 2D with FTCS + * approach and constant boundary condition + * */ #include // #include // #define EASY_PROFILER_ENABLE ::profiler::setEnabled(true); - int main(int argc, char *argv[]) { - // EASY_PROFILER_ENABLE; - // profiler::startListen(); - // ************** - // **** GRID **** - // ************** - // profiler::startListen(); - // create a grid with a 20 x 20 field - int row = 20; - int col = 20; - Grid grid = Grid(row,col); + // EASY_PROFILER_ENABLE; + // profiler::startListen(); + // ************** + // **** GRID **** + // ************** + // profiler::startListen(); + // create a grid with a 20 x 20 field + int row = 20; + int col = 20; + Grid grid = Grid(row, col); - // (optional) set the domain, e.g.: - // grid.setDomain(20, 20); + // (optional) set the domain, e.g.: + // grid.setDomain(20, 20); - // (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,0); - concentrations(0,0) = 1999; - grid.setConcentrations(concentrations); - + // (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, 0); + concentrations(0, 0) = 1999; + grid.setConcentrations(concentrations); - // (optional) set alphas of the grid, e.g.: - // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value - // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value - // grid.setAlpha(alphax, alphay); + // (optional) set alphas of the grid, e.g.: + // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value + // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value + // grid.setAlpha(alphax, alphay); + // ****************** + // **** BOUNDARY **** + // ****************** - // ****************** - // **** BOUNDARY **** - // ****************** + // 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); - // 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); + // (optional) set boundary condition values for one side, e.g.: + // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value + // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values + // VectorXd bc_zero_values = VectorXd::Constant(20,0); + // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); + // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); + // VectorXd bc_front_values = VectorXd::Constant(20,2000); + // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); + // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); + // ************************ + // **** SIMULATION ENV **** + // ************************ - // (optional) set boundary condition values for one side, e.g.: - // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value - // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values - // VectorXd bc_zero_values = VectorXd::Constant(20,0); - // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values); - // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values); - // VectorXd bc_front_values = VectorXd::Constant(20,2000); - // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values); - // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values); + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + // set the timestep of the simulation + simulation.setTimestep(0.1); // timestep - // ************************ - // **** SIMULATION ENV **** - // ************************ + // set the number of iterations + simulation.setIterations(10000); - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - // set the timestep of the simulation - simulation.setTimestep(0.1); // timestep + // **** RUN SIMULATION **** - // set the number of iterations - simulation.setIterations(10000); + // run the simulation - // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); - - - // **** RUN SIMULATION **** - - // run the simulation - - // EASY_BLOCK("SIMULATION") - simulation.run(); - // EASY_END_BLOCK; - // profiler::dumpBlocksToFile("test_profile.prof"); - // profiler::stopListen(); + // EASY_BLOCK("SIMULATION") + simulation.run(); + // EASY_END_BLOCK; + // profiler::dumpBlocksToFile("test_profile.prof"); + // profiler::stopListen(); } \ No newline at end of file diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index 80d4112..9189555 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -1,77 +1,78 @@ /** * @file FTCS_2D_proto_example.cpp * @author Hannes Signer, Philipp Ungrund - * @brief Creates a prototypical standard TUG simulation in 2D with FTCS approach - * and constant boundary condition - * + * @brief Creates a prototypical standard TUG simulation in 2D with FTCS + * approach and constant boundary condition + * */ #include int main(int argc, char *argv[]) { - - // ************** - // **** GRID **** - // ************** - // create a grid with a 20 x 20 field - int row = 64; - int col = 64; - int n2 = row/2-1; - Grid grid = Grid(row,col); + // ************** + // **** GRID **** + // ************** - // (optional) set the domain, e.g.: - // grid.setDomain(20, 20); + // create a grid with a 20 x 20 field + int row = 64; + int col = 64; + int n2 = row / 2 - 1; + Grid grid = Grid(row, col); - // (optional) set the concentrations, e.g.: - // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value - MatrixXd concentrations = MatrixXd::Constant(row,col,0); - concentrations(n2,n2) = 1; - concentrations(n2,n2+1) = 1; - concentrations(n2+1,n2) = 1; - concentrations(n2+1,n2+1) = 1; - grid.setConcentrations(concentrations); + // (optional) set the domain, e.g.: + // grid.setDomain(20, 20); - // (optional) set alphas of the grid, e.g.: - MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value - MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value - grid.setAlpha(alphax, alphay); + // (optional) set the concentrations, e.g.: + // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // + // #row,#col,value + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(n2, n2) = 1; + concentrations(n2, n2 + 1) = 1; + concentrations(n2 + 1, n2) = 1; + concentrations(n2 + 1, n2 + 1) = 1; + grid.setConcentrations(concentrations); + // (optional) set alphas of the grid, e.g.: + MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value + MatrixXd alphay = MatrixXd::Constant(row, col, 1E-6); // row,col,value + grid.setAlpha(alphax, alphay); - // ****************** - // **** BOUNDARY **** - // ****************** + // ****************** + // **** BOUNDARY **** + // ****************** - // create a boundary with constant values - Boundary bc = Boundary(grid); + // create a boundary with constant values + Boundary bc = Boundary(grid); - // (optional) set boundary condition values for one side, e.g.: - bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); // side,values - 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.: + bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); // side,values + bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + bc.setBoundarySideConstant(BC_SIDE_TOP, 0); + bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); + // ************************ + // **** SIMULATION ENV **** + // ************************ - // ************************ - // **** SIMULATION ENV **** - // ************************ + // set up a simulation environment + Simulation simulation = + Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach - // set up a simulation environment - Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + // (optional) set the timestep of the simulation + simulation.setTimestep(1000); // timestep - // (optional) set the timestep of the simulation - simulation.setTimestep(1000); // timestep + // (optional) set the number of iterations + simulation.setIterations(5); - // (optional) set the number of iterations - simulation.setIterations(5); + // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_OFF); - // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE] - simulation.setOutputCSV(CSV_OUTPUT_OFF); - - // **** RUN SIMULATION **** + // **** RUN SIMULATION **** - // run the simulation - simulation.run(); + // run the simulation + simulation.run(); - return 0; + return 0; } diff --git a/examples/profiling_openmp.cpp b/examples/profiling_openmp.cpp index 78c915b..92b5b12 100644 --- a/examples/profiling_openmp.cpp +++ b/examples/profiling_openmp.cpp @@ -1,67 +1,66 @@ +#include +#include +#include #include #include -#include -#include -#include int main(int argc, char *argv[]) { - int n[] = {2000}; - int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; - int iterations[1] = {1}; - int repetition = 10; + int n[] = {2000}; + int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; + int iterations[1] = {1}; + int repetition = 10; - - for(int l=0; l(end - begin); - myfile << milliseconds.count() << endl; - } + auto begin = std::chrono::high_resolution_clock::now(); + sim.run(); + auto end = std::chrono::high_resolution_clock::now(); + auto milliseconds = + std::chrono::duration_cast(end - + begin); + myfile << milliseconds.count() << endl; } - cout << endl; - myfile << endl; - + } + cout << endl; + myfile << endl; } myfile.close(); - } + } } \ No newline at end of file diff --git a/examples/profiling_speedup.cpp b/examples/profiling_speedup.cpp index 78c915b..92b5b12 100644 --- a/examples/profiling_speedup.cpp +++ b/examples/profiling_speedup.cpp @@ -1,67 +1,66 @@ +#include +#include +#include #include #include -#include -#include -#include int main(int argc, char *argv[]) { - int n[] = {2000}; - int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; - int iterations[1] = {1}; - int repetition = 10; + int n[] = {2000}; + int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; + int iterations[1] = {1}; + int repetition = 10; - - for(int l=0; l(end - begin); - myfile << milliseconds.count() << endl; - } + auto begin = std::chrono::high_resolution_clock::now(); + sim.run(); + auto end = std::chrono::high_resolution_clock::now(); + auto milliseconds = + std::chrono::duration_cast(end - + begin); + myfile << milliseconds.count() << endl; } - cout << endl; - myfile << endl; - + } + cout << endl; + myfile << endl; } myfile.close(); - } + } } \ No newline at end of file diff --git a/examples/reference-FTCS_2D_closed.cpp b/examples/reference-FTCS_2D_closed.cpp index 7b51e67..09dad0d 100644 --- a/examples/reference-FTCS_2D_closed.cpp +++ b/examples/reference-FTCS_2D_closed.cpp @@ -1,55 +1,51 @@ -#include #include "Eigen/Core" #include +#include using namespace std; int main(int argc, char *argv[]) { - int row = 50; - int col = 50; - int domain_row = 10; - int domain_col = 10; + int row = 50; + int col = 50; + int domain_row = 10; + int domain_col = 10; + // Grid + Grid grid = Grid(row, col); + grid.setDomain(domain_row, domain_col); - // Grid - Grid grid = Grid(row, col); - grid.setDomain(domain_row, domain_col); + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(5, 5) = 1; + grid.setConcentrations(concentrations); - MatrixXd concentrations = MatrixXd::Constant(row, col, 0); - concentrations(5,5) = 1; - grid.setConcentrations(concentrations); - - MatrixXd alpha = MatrixXd::Constant(row, col, 1); - for (int i = 0; i < 5; i++) { - for (int j = 0; j < 6; j++) { - alpha(i, j) = 0.01; - } + MatrixXd alpha = MatrixXd::Constant(row, col, 1); + for (int i = 0; i < 5; i++) { + for (int j = 0; j < 6; j++) { + alpha(i, j) = 0.01; } - for (int i = 0; i < 5; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.001; - } + } + for (int i = 0; i < 5; i++) { + for (int j = 6; j < 11; j++) { + alpha(i, j) = 0.001; } - for (int i = 5; i < 11; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.1; - } + } + for (int i = 5; i < 11; i++) { + for (int j = 6; j < 11; j++) { + alpha(i, j) = 0.1; } - grid.setAlpha(alpha, alpha); + } + grid.setAlpha(alpha, alpha); + // Boundary + Boundary bc = Boundary(grid); - // Boundary - Boundary bc = Boundary(grid); + // Simulation + Simulation sim = Simulation(grid, bc, FTCS_APPROACH); + sim.setTimestep(0.001); + sim.setIterations(10000); + sim.setOutputCSV(CSV_OUTPUT_OFF); + sim.setOutputConsole(CONSOLE_OUTPUT_OFF); - - // Simulation - Simulation sim = Simulation(grid, bc, FTCS_APPROACH); - sim.setTimestep(0.001); - sim.setIterations(10000); - sim.setOutputCSV(CSV_OUTPUT_OFF); - sim.setOutputConsole(CONSOLE_OUTPUT_OFF); - - - // RUN - sim.run(); + // RUN + sim.run(); } \ No newline at end of file diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index 319503a..3b900e3 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -1,37 +1,29 @@ /** * @file Boundary.hpp - * @brief API of Boundary class, that holds all information for each boundary condition - * at the edges of the diffusion grid. - * + * @brief API of Boundary class, that holds all information for each boundary + * condition at the edges of the diffusion grid. + * */ #ifndef BOUNDARY_H_ #define BOUNDARY_H_ -#include #include "Grid.hpp" +#include using namespace std; using namespace Eigen; /** - * @brief Enum defining the two implemented boundary conditions. - * + * @brief Enum defining the two implemented boundary conditions. + * */ -enum BC_TYPE { - BC_TYPE_CLOSED, - BC_TYPE_CONSTANT -}; +enum BC_TYPE { BC_TYPE_CLOSED, BC_TYPE_CONSTANT }; /** - * @brief Enum defining all 4 possible sides to a 1D and 2D grid. - * + * @brief Enum defining all 4 possible sides to a 1D and 2D grid. + * */ -enum BC_SIDE { - BC_SIDE_LEFT, - BC_SIDE_RIGHT, - BC_SIDE_TOP, - BC_SIDE_BOTTOM -}; +enum BC_SIDE { BC_SIDE_LEFT, BC_SIDE_RIGHT, BC_SIDE_TOP, BC_SIDE_BOTTOM }; /** * This class defines the boundary conditions of individual boundary elements. @@ -39,177 +31,184 @@ enum BC_SIDE { * The class serves as an auxiliary class for structuring the Boundary class. */ class BoundaryElement { - public: - - /** - * @brief Construct a new Boundary Element object for the closed case. - * The boundary type is here automatically set to the type - * BC_TYPE_CLOSED, where the value takes -1 and does not hold any - * physical meaning. - */ - BoundaryElement(); +public: + /** + * @brief Construct a new Boundary Element object for the closed case. + * The boundary type is here automatically set to the type + * BC_TYPE_CLOSED, where the value takes -1 and does not hold any + * physical meaning. + */ + BoundaryElement(); - /** - * @brief Construct a new Boundary Element object for the constant case. - * The boundary type is automatically set to the type - * BC_TYPE_CONSTANT. - * - * @param value Value of the constant concentration to be assumed at the - * corresponding boundary element. - */ - BoundaryElement(double value); + /** + * @brief Construct a new Boundary Element object for the constant case. + * The boundary type is automatically set to the type + * BC_TYPE_CONSTANT. + * + * @param value Value of the constant concentration to be assumed at the + * corresponding boundary element. + */ + BoundaryElement(double value); - /** - * @brief Allows changing the boundary type of a corresponding - * BoundaryElement object. - * - * @param type Type of boundary condition. Either BC_TYPE_CONSTANT or - BC_TYPE_CLOSED. - */ - void setType(BC_TYPE type); - - /** - * @brief Sets the value of a boundary condition for the constant case. - * - * @param value Concentration to be considered constant for the - * corresponding boundary element. - */ - void setValue(double value); + /** + * @brief Allows changing the boundary type of a corresponding + * BoundaryElement object. + * + * @param type Type of boundary condition. Either BC_TYPE_CONSTANT or + BC_TYPE_CLOSED. + */ + void setType(BC_TYPE type); - /** - * @brief Return the type of the boundary condition, i.e. whether the - * boundary is considered closed or constant. - * - * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or - BC_TYPE_CONSTANT. - */ - BC_TYPE getType(); + /** + * @brief Sets the value of a boundary condition for the constant case. + * + * @param value Concentration to be considered constant for the + * corresponding boundary element. + */ + void setValue(double value); - /** - * @brief Return the concentration value for the constant boundary condition. - * - * @return double Value of the concentration. - */ - double getValue(); + /** + * @brief Return the type of the boundary condition, i.e. whether the + * boundary is considered closed or constant. + * + * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or + BC_TYPE_CONSTANT. + */ + BC_TYPE getType(); - private: - BC_TYPE type; - double value; + /** + * @brief Return the concentration value for the constant boundary condition. + * + * @return double Value of the concentration. + */ + double getValue(); + +private: + BC_TYPE type; + double value; }; - /** * This class implements the functionality and management of the boundary * conditions in the grid to be simulated. */ class Boundary { - public: - /** - * @brief Creates a boundary object based on the passed grid object and - * initializes the boundaries as closed. - * - * @param grid Grid object on the basis of which the simulation takes place - * and from which the dimensions (in 2D case) are taken. - */ - Boundary(Grid grid); +public: + /** + * @brief Creates a boundary object based on the passed grid object and + * initializes the boundaries as closed. + * + * @param grid Grid object on the basis of which the simulation takes place + * and from which the dimensions (in 2D case) are taken. + */ + Boundary(Grid grid); - /** - * @brief Sets all elements of the specified boundary side to the boundary - * condition closed. - * - * @param side Side to be set to closed, e.g. BC_SIDE_LEFT. - */ - void setBoundarySideClosed(BC_SIDE side); + /** + * @brief Sets all elements of the specified boundary side to the boundary + * condition closed. + * + * @param side Side to be set to closed, e.g. BC_SIDE_LEFT. + */ + void setBoundarySideClosed(BC_SIDE side); - /** - * @brief Sets all elements of the specified boundary side to the boundary - * condition constant. Thereby the concentration values of the - * boundaries are set to the passed value. - * - * @param side Side to be set to constant, e.g. BC_SIDE_LEFT. - * @param value Concentration to be set for all elements of the specified page. - */ - void setBoundarySideConstant(BC_SIDE side, double value); + /** + * @brief Sets all elements of the specified boundary side to the boundary + * condition constant. Thereby the concentration values of the + * boundaries are set to the passed value. + * + * @param side Side to be set to constant, e.g. BC_SIDE_LEFT. + * @param value Concentration to be set for all elements of the specified + * page. + */ + void setBoundarySideConstant(BC_SIDE side, double value); - /** - * @brief Specifically sets the boundary element of the specified side - * defined by the index to the boundary condition closed. - * - * @param side Side in which an element is to be defined as closed. - * @param index Index of the boundary element on the corresponding - * boundary side. Must index an element of the corresponding side. - */ - void setBoundaryElementClosed(BC_SIDE side, int index); + /** + * @brief Specifically sets the boundary element of the specified side + * defined by the index to the boundary condition closed. + * + * @param side Side in which an element is to be defined as closed. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + * side. + */ + void setBoundaryElementClosed(BC_SIDE side, int index); - /** - * @brief Specifically sets the boundary element of the specified side - * defined by the index to the boundary condition constant with the - given concentration value. - * - * @param side Side in which an element is to be defined as constant. - * @param index Index of the boundary element on the corresponding - * boundary side. Must index an element of the corresponding side. - * @param value Concentration value to which the boundary element should be set. - */ - void setBoundaryElementConstant(BC_SIDE side, int index, double value); + /** + * @brief Specifically sets the boundary element of the specified side + * defined by the index to the boundary condition constant with the + given concentration value. + * + * @param side Side in which an element is to be defined as constant. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + side. + * @param value Concentration value to which the boundary element should be + set. + */ + void setBoundaryElementConstant(BC_SIDE side, int index, double value); - /** - * @brief Returns the boundary condition of a specified side as a vector - * of BoundarsElement objects. - * - * @param side Boundary side from which the boundary conditions are to be returned. - * @return vector Contains the boundary conditions as BoundaryElement objects. - */ - const vector getBoundarySide(BC_SIDE side); + /** + * @brief Returns the boundary condition of a specified side as a vector + * of BoundarsElement objects. + * + * @param side Boundary side from which the boundary conditions are to be + * returned. + * @return vector Contains the boundary conditions as + * BoundaryElement objects. + */ + const vector getBoundarySide(BC_SIDE side); - /** - * @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); + /** + * @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); - /** - * @brief Returns the boundary condition of a specified element on a given side. - * - * @param side Boundary side in which the boundary condition is located. - * @param index Index of the boundary element on the corresponding - * boundary side. Must index an element of the corresponding side. - * @return BoundaryElement Boundary condition as a BoundaryElement object. - */ - BoundaryElement getBoundaryElement(BC_SIDE side, int index); + /** + * @brief Returns the boundary condition of a specified element on a given + * side. + * + * @param side Boundary side in which the boundary condition is located. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + * side. + * @return BoundaryElement Boundary condition as a BoundaryElement object. + */ + BoundaryElement getBoundaryElement(BC_SIDE side, int index); - /** - * @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED or - BC_TYPE_CONSTANT. - * - * @param side Boundary side in which the boundary condition type is located. - * @param index Index of the boundary element on the corresponding - * boundary side. Must index an element of the corresponding side. - * @return BC_TYPE Boundary Type of the corresponding boundary condition. - */ - BC_TYPE getBoundaryElementType(BC_SIDE side, int index); + /** + * @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED + or BC_TYPE_CONSTANT. + * + * @param side Boundary side in which the boundary condition type is located. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + side. + * @return BC_TYPE Boundary Type of the corresponding boundary condition. + */ + BC_TYPE getBoundaryElementType(BC_SIDE side, int index); - /** - * @brief Returns the concentration value of a corresponding - * BoundaryElement object if it is a constant boundary condition. - * - * @param side Boundary side in which the boundary condition value is - * located. - * @param index Index of the boundary element on the corresponding - * boundary side. Must index an element of the corresponding - * side. - * @return double Concentration of the corresponding BoundaryElement object. - */ - double getBoundaryElementValue(BC_SIDE side, int index); + /** + * @brief Returns the concentration value of a corresponding + * BoundaryElement object if it is a constant boundary condition. + * + * @param side Boundary side in which the boundary condition value is + * located. + * @param index Index of the boundary element on the corresponding + * boundary side. Must index an element of the corresponding + * side. + * @return double Concentration of the corresponding BoundaryElement object. + */ + double getBoundaryElementValue(BC_SIDE side, int index); - private: - Grid grid; // Boundary is directly dependent on the dimensions of a predefined - - vector> boundaries; // Vector with Boundary Element information +private: + Grid grid; // Boundary is directly dependent on the dimensions of a predefined + + vector> + boundaries; // Vector with Boundary Element information }; #endif - diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index c7c343e..09ae737 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -1,8 +1,8 @@ /** * @file Grid.hpp - * @brief API of Grid class, that holds a matrix with concenctrations and a - * respective matrix/matrices of alpha coefficients. - * + * @brief API of Grid class, that holds a matrix with concenctrations and a + * respective matrix/matrices of alpha coefficients. + * */ #include @@ -11,165 +11,166 @@ using namespace Eigen; class Grid { - public: +public: + /** + * @brief Constructs a new 1D-Grid object of a given length, which holds a + * matrix with concentrations and a respective matrix of alpha coefficients. + * The domain length is per default the same as the length. The + * concentrations are all 20 by default and the alpha coefficients are 1. + * + * @param length Length of the 1D-Grid. Must be greater than 3. + */ + Grid(int length); - /** - * @brief Constructs a new 1D-Grid object of a given length, which holds a matrix - * with concentrations and a respective matrix of alpha coefficients. - * The domain length is per default the same as the length. The concentrations - * are all 20 by default and the alpha coefficients are 1. - * - * @param length Length of the 1D-Grid. Must be greater than 3. - */ - Grid(int length); + /** + * @brief Constructs a new 2D-Grid object of given dimensions, which holds a + * matrix with concentrations and the respective matrices of alpha coefficient + * for each direction. The domain in x- and y-direction is per default equal + * to the col length and row length, respectively. The concentrations are all + * 20 by default across the entire grid and the alpha coefficients 1 in both + * directions. + * + * @param row Length of the 2D-Grid in y-direction. Must be greater than 3. + * @param col Length of the 2D-Grid in x-direction. Must be greater than 3. + */ + Grid(int row, int col); - /** - * @brief Constructs a new 2D-Grid object of given dimensions, which holds a matrix - * with concentrations and the respective matrices of alpha coefficient for - * each direction. The domain in x- and y-direction is per default equal to - * the col length and row length, respectively. - * The concentrations are all 20 by default across the entire grid and the - * alpha coefficients 1 in both directions. - * - * @param row Length of the 2D-Grid in y-direction. Must be greater than 3. - * @param col Length of the 2D-Grid in x-direction. Must be greater than 3. - */ - Grid(int row, int col); + /** + * @brief Sets the concentrations matrix for a 1D or 2D-Grid. + * + * @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix + * must have correct dimensions as defined in row and col. (Or length, in 1D + * case). + */ + void setConcentrations(MatrixXd concentrations); - /** - * @brief Sets the concentrations matrix for a 1D or 2D-Grid. - * - * @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix must - * have correct dimensions as defined in row and col. (Or length, - * in 1D case). - */ - void setConcentrations(MatrixXd concentrations); + /** + * @brief Gets the concentrations matrix for a Grid. + * + * @return MatrixXd An Eigen3 matrix holding the concentrations and having the + * same dimensions as the grid. + */ + const MatrixXd getConcentrations(); - /** - * @brief Gets the concentrations matrix for a Grid. - * - * @return MatrixXd An Eigen3 matrix holding the concentrations and having the - * same dimensions as the grid. - */ - const MatrixXd getConcentrations(); + /** + * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one + * dimensional. + * + * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients. + * Matrix columns must have same size as length of grid. + */ + void setAlpha(MatrixXd alpha); - /** - * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one dimensional. - * - * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients. Matrix - * columns must have same size as length of grid. - */ - void setAlpha(MatrixXd alpha); + /** + * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two + * dimensional. + * + * @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in + * x-direction. Matrix must be of same size as the grid. + * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in + * y-direction. Matrix must be of same size as the grid. + */ + void setAlpha(MatrixXd alphaX, MatrixXd alphaY); - /** - * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two dimensional. - * - * @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in x-direction. - * Matrix must be of same size as the grid. - * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in y-direction. - * Matrix must be of same size as the grid. - */ - void setAlpha(MatrixXd alphaX, MatrixXd alphaY); + /** + * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one + * dimensional. + * + * @return MatrixXd A matrix with 1 row holding the alpha coefficients. + */ + const MatrixXd getAlpha(); - /** - * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one dimensional. - * - * @return MatrixXd A matrix with 1 row holding the alpha coefficients. - */ - const MatrixXd getAlpha(); + /** + * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. + * Grid must be two dimensional. + * + * @return MatrixXd A matrix holding the alpha coefficients in x-direction. + */ + const MatrixXd getAlphaX(); - /** - * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. Grid must be - * two dimensional. - * - * @return MatrixXd A matrix holding the alpha coefficients in x-direction. - */ - const MatrixXd getAlphaX(); + /** + * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. + * Grid must be two dimensional. + * + * @return MatrixXd A matrix holding the alpha coefficients in y-direction. + */ + const MatrixXd getAlphaY(); - /** - * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. Grid must be - * two dimensional. - * - * @return MatrixXd A matrix holding the alpha coefficients in y-direction. - */ - const MatrixXd getAlphaY(); + /** + * @brief Gets the dimensions of the grid. + * + * @return int Dimensions, either 1 or 2. + */ + int getDim(); - /** - * @brief Gets the dimensions of the grid. - * - * @return int Dimensions, either 1 or 2. - */ - int getDim(); + /** + * @brief Gets length of 1D grid. Must be one dimensional grid. + * + * @return int Length of 1D grid. + */ + int getLength(); - /** - * @brief Gets length of 1D grid. Must be one dimensional grid. - * - * @return int Length of 1D grid. - */ - int getLength(); + /** + * @brief Gets the number of rows of the grid. + * + * @return int Number of rows. + */ + int getRow(); - /** - * @brief Gets the number of rows of the grid. - * - * @return int Number of rows. - */ - int getRow(); + /** + * @brief Gets the number of columns of the grid. + * + * @return int Number of columns. + */ + int getCol(); - /** - * @brief Gets the number of columns of the grid. - * - * @return int Number of columns. - */ - int getCol(); + /** + * @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional. + * + * @param domainLength A double value of the domain length. Must be positive. + */ + void setDomain(double domainLength); - /** - * @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional. - * - * @param domainLength A double value of the domain length. Must be positive. - */ - void setDomain(double domainLength); + /** + * @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional. + * + * @param domainRow A double value of the domain size in y-direction. Must be + * positive. + * @param domainCol A double value of the domain size in x-direction. Must be + * positive. + */ + void setDomain(double domainRow, double domainCol); - /** - * @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional. - * - * @param domainRow A double value of the domain size in y-direction. Must be positive. - * @param domainCol A double value of the domain size in x-direction. Must be positive. - */ - void setDomain(double domainRow,double 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 for 1D-Grid. Grid must be one dimensional. - * - * @return double Delta value. - */ - double getDelta(); + /** + * @brief Gets the delta value in x-direction. + * + * @return double Delta value in x-direction. + */ + double getDeltaCol(); - /** - * @brief Gets the delta value in x-direction. - * - * @return double Delta value in x-direction. - */ - double getDeltaCol(); - - /** - * @brief Gets the delta value in y-direction. Must be two dimensional grid. - * - * @return double Delta value in y-direction. - */ - double getDeltaRow(); - - - private: - - int col; // number of grid columns - int row; // number of grid rows - int dim; // 1D or 2D - double domainCol; // number of domain columns - double domainRow; // number of domain rows - double deltaCol; // delta in x-direction (between columns) - double deltaRow; // delta in y-direction (between rows) - MatrixXd concentrations; // Matrix holding grid concentrations - MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction - MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction + /** + * @brief Gets the delta value in y-direction. Must be two dimensional grid. + * + * @return double Delta value in y-direction. + */ + double getDeltaRow(); +private: + int col; // number of grid columns + int row; // number of grid rows + int dim; // 1D or 2D + double domainCol; // number of domain columns + double domainRow; // number of domain rows + double deltaCol; // delta in x-direction (between columns) + double deltaRow; // delta in y-direction (between rows) + MatrixXd concentrations; // Matrix holding grid concentrations + MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction + MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction }; diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 641167b..f74ada6 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -1,8 +1,8 @@ /** * @file Simulation.hpp - * @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. + * @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" @@ -11,52 +11,54 @@ using namespace std; /** - * @brief Enum defining the two implemented solution approaches. - * + * @brief Enum defining the two implemented solution approaches. + * */ enum APPROACH { - FTCS_APPROACH, // Forward Time-Centered Space - BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver - CRANK_NICOLSON_APPROACH + FTCS_APPROACH, // Forward Time-Centered Space + BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver + CRANK_NICOLSON_APPROACH }; /** * @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 + EIGEN_LU_SOLVER, // EigenLU solver + THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for + // tridiagonal matrices }; /** * @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 // csv output like VERBOSE but additional boundary conditions at beginning + 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 // 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 + 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. - * + * @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_OFF, // do not print any time measures + TIME_MEASURE_ON // print time measure after last iteration }; /** @@ -66,148 +68,155 @@ enum TIME_MEASURE { * */ class Simulation { - public: - /** - * @brief Set up a simulation environment. The timestep and number of iterations - * must be set. For the BTCS approach, the Thomas algorithm is used as - * the default linear equation solver as this is faster for tridiagonal - * matrices. CSV output, console output and time measure are off by default. - * Also, the number of cores is set to the maximum number of cores -1 by default. - * - * @param grid Valid grid object - * @param bc Valid boundary condition object - * @param approach Approach to solving the problem. Either FTCS or BTCS. - */ - Simulation(Grid &grid, Boundary &bc, APPROACH approach); +public: + /** + * @brief Set up a simulation environment. The timestep and number of + * iterations must be set. For the BTCS approach, the Thomas algorithm is used + * as the default linear equation solver as this is faster for tridiagonal + * matrices. CSV output, console output and time measure are off by + * default. Also, the number of cores is set to the maximum number of cores -1 + * by default. + * + * @param grid Valid grid object + * @param bc Valid boundary condition object + * @param approach Approach to solving the problem. Either FTCS or BTCS. + */ + Simulation(Grid &grid, Boundary &bc, APPROACH approach); - /** - * @brief Set the option to output the results to a CSV file. Off by default. - * - * - * @param csv_output Valid output option. The following options can be set - * here: - * - 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 - */ - void setOutputCSV(CSV_OUTPUT csv_output); + /** + * @brief Set the option to output the results to a CSV file. Off by default. + * + * + * @param csv_output Valid output option. The following options can be set + * here: + * - 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 + */ + void setOutputCSV(CSV_OUTPUT csv_output); - /** - * @brief Set the options for outputting information to the console. Off by default. - * - * @param console_output Valid output option. The following options can be set - * here: - * - 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 - */ - void setOutputConsole(CONSOLE_OUTPUT console_output); + /** + * @brief Set the options for outputting information to the console. Off by + * default. + * + * @param console_output Valid output option. The following options can be set + * here: + * - 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 + */ + void setOutputConsole(CONSOLE_OUTPUT console_output); - /** - * @brief Set the Time Measure option. Off by default. - * - * @param time_measure The following options are allowed: - * - TIME_MEASURE_OFF: Time of simulation is not printed to console - * - TIME_MEASURE_ON: Time of simulation run is printed to console - */ - void setTimeMeasure(TIME_MEASURE time_measure); + /** + * @brief Set the Time Measure option. Off by default. + * + * @param time_measure The following options are allowed: + * - TIME_MEASURE_OFF: Time of simulation is not printed + * to console + * - TIME_MEASURE_ON: Time of simulation run is printed to + * console + */ + void setTimeMeasure(TIME_MEASURE time_measure); - /** - * @brief Setting the time step for each iteration step. Time step must be - * greater than zero. Setting the timestep is required. - * - * @param timestep Valid timestep greater than zero. - */ - void setTimestep(double timestep); + /** + * @brief Setting the time step for each iteration step. Time step must be + * greater than zero. Setting the timestep is required. + * + * @param timestep Valid timestep greater than zero. + */ + void setTimestep(double timestep); - /** - * @brief Currently set time step is returned. - * - * @return double timestep - */ - double getTimestep(); + /** + * @brief Currently set time step is returned. + * + * @return double timestep + */ + double getTimestep(); - /** - * @brief Set the desired iterations to be calculated. A value greater - * than zero must be specified here. Setting iterations is required. - * - * @param iterations Number of iterations to be simulated. - */ - void setIterations(int iterations); + /** + * @brief Set the desired iterations to be calculated. A value greater + * than zero must be specified here. Setting iterations is required. + * + * @param iterations Number of iterations to be simulated. + */ + 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 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 Set the number of desired openMP Threads. - * - * @param num_threads Number of desired threads. Must have a value between - * 1 and the maximum available number of processors. The maximum number of - * processors is set as the default case during Simulation construction. - */ - void setNumberThreads(int num_threads); + /** + * @brief Set the number of desired openMP Threads. + * + * @param num_threads Number of desired threads. Must have a value between + * 1 and the maximum available number of processors. The + * maximum number of processors is set as the default case during Simulation + * construction. + */ + void setNumberThreads(int num_threads); - /** - * @brief Return the currently set iterations to be calculated. - * - * @return int Number of iterations. - */ - int getIterations(); + /** + * @brief Return the currently set iterations to be calculated. + * + * @return int Number of iterations. + */ + int getIterations(); - /** - * @brief Outputs the current concentrations of the grid on the console. - * - */ - void printConcentrationsConsole(); + /** + * @brief Outputs the current concentrations of the grid on the console. + * + */ + void printConcentrationsConsole(); - /** - * @brief Creates a CSV file with a name containing the current simulation - * parameters. If the data name already exists, an additional counter is - * appended to the name. The name of the file is built up as follows: - * + + + +.csv - * - * @return string Filename with configured simulation parameters. - */ - string createCSVfile(); + /** + * @brief Creates a CSV file with a name containing the current simulation + * parameters. If the data name already exists, an additional counter + * is appended to the name. The name of the file is built up as follows: + * + + + +.csv + * + * @return string Filename with configured simulation parameters. + */ + string createCSVfile(); - /** - * @brief Writes the currently calculated concentration values of the grid - * into the CSV file with the passed filename. - * - * @param filename Name of the file to which the concentration values are - * to be written. - */ - void printConcentrationsCSV(string filename); + /** + * @brief Writes the currently calculated concentration values of the grid + * into the CSV file with the passed filename. + * + * @param filename Name of the file to which the concentration values are + * to be written. + */ + void printConcentrationsCSV(string filename); - /** - * @brief Method starts the simulation process with the previously set - * parameters. - */ - void run(); + /** + * @brief Method starts the simulation process with the previously set + * parameters. + */ + void run(); - private: - - double timestep; - int iterations; - int innerIterations; - int numThreads; - CSV_OUTPUT csv_output; - CONSOLE_OUTPUT console_output; - TIME_MEASURE time_measure; - - Grid &grid; - Boundary &bc; - APPROACH approach; - SOLVER solver; +private: + double timestep; + int iterations; + int innerIterations; + int numThreads; + CSV_OUTPUT csv_output; + CONSOLE_OUTPUT console_output; + TIME_MEASURE time_measure; + Grid &grid; + Boundary &bc; + APPROACH approach; + SOLVER solver; }; diff --git a/src/BTCSv2.cpp b/src/BTCSv2.cpp index 49e4d1c..2ff52b8 100644 --- a/src/BTCSv2.cpp +++ b/src/BTCSv2.cpp @@ -1,437 +1,459 @@ /** * @file BTCSv2.cpp - * @brief Implementation of heterogenous BTCS (backward time-centered space) solution - * of diffusion equation in 1D and 2D space. Internally the alternating-direction - * implicit (ADI) method is used. Version 2, because Version 1 was an - * implementation for the homogeneous BTCS solution. - * + * @brief Implementation of heterogenous BTCS (backward time-centered space) + * solution of diffusion equation in 1D and 2D space. Internally the + * alternating-direction implicit (ADI) method is used. Version 2, because + * Version 1 was an implementation for the homogeneous BTCS solution. + * */ #include "FTCS.cpp" -#include #include +#include #define NUM_THREADS_BTCS 10 using namespace Eigen; - // calculates coefficient for left boundary in constant case -static tuple calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) { - double centerCoeff; - double rightCoeff; +static tuple +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)); + 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}; + return {centerCoeff, rightCoeff}; } - // calculates coefficient for left boundary in closed case -static tuple calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) { - double centerCoeff; - double rightCoeff; +static tuple +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)); + centerCoeff = + 1 + sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)); + rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)); - return {centerCoeff, rightCoeff}; + return {centerCoeff, rightCoeff}; } - // calculates coefficient for right boundary in constant case -static tuple calcRightBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, int n, double sx) { - double leftCoeff; - double centerCoeff; +static tuple 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)); + 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}; + return {leftCoeff, centerCoeff}; } - // calculates coefficient for right boundary in closed case -static tuple calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) { - double leftCoeff; - double centerCoeff; +static tuple +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)); + leftCoeff = + -sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); + centerCoeff = + 1 + sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); - return {leftCoeff, centerCoeff}; + return {leftCoeff, centerCoeff}; } - // creates coefficient matrix for next time step from alphas in x-direction -static SparseMatrix createCoeffMatrix(MatrixXd &alpha, vector &bcLeft, vector &bcRight, int numCols, int rowIndex, double sx) { +static SparseMatrix createCoeffMatrix(MatrixXd &alpha, + vector &bcLeft, + vector &bcRight, + int numCols, int rowIndex, + double sx) { - // square matrix of column^2 dimension for the coefficients - SparseMatrix cm(numCols, numCols); - cm.reserve(VectorXi::Constant(numCols, 3)); + // square matrix of column^2 dimension for the coefficients + SparseMatrix cm(numCols, numCols); + cm.reserve(VectorXi::Constant(numCols, 3)); - // left column - 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!"); - } + // left column + 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; - for (int i = 1; i < n; i++) { - cm.insert(i,i-1) = -sx * calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i)); - cm.insert(i,i) = 1 + sx * ( - calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex,i+1)) - + calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i)) - ) - ; - cm.insert(i,i+1) = -sx * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex,i+1)); - } + // inner columns + int n = numCols - 1; + for (int i = 1; i < n; i++) { + cm.insert(i, i - 1) = + -sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i)); + cm.insert(i, i) = + 1 + + sx * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)) + + calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i))); + cm.insert(i, i + 1) = + -sx * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)); + } - // right column - 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!"); - } + // right column + 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(); // important for Eigen solver + cm.makeCompressed(); // important for Eigen solver - return cm; + return cm; } - // calculates explicity concentration at top boundary in constant case -static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrations, - MatrixXd &alpha, vector &bcTop, int rowIndex, int i, double sy) { - double c; +static double calcExplicitConcentrationsTopBoundaryConstant( + MatrixXd &concentrations, MatrixXd &alpha, vector &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(); + 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; + 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; +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); + 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; + return c; } - // calculates explicit concentration at bottom boundary in constant case -static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concentrations, - MatrixXd &alpha, vector &bcBottom, int rowIndex, int i, double sy) { - double c; +static double calcExplicitConcentrationsBottomBoundaryConstant( + MatrixXd &concentrations, MatrixXd &alpha, + vector &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); + 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; + 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; +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); + 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; + 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, + vector &bcLeft, vector &bcRight, + vector &bcTop, vector &bcBottom, + int length, int rowIndex, double sx, double sy) { -// creates a solution vector for next time step from the current state of concentrations -static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, - vector &bcLeft, vector &bcRight, - vector &bcTop, vector &bcBottom, - int length, int rowIndex, double sx, double sy) { + VectorXd sv(length); + int numRows = concentrations.rows(); + BC_TYPE type; - VectorXd sv(length); - int numRows = concentrations.rows(); - BC_TYPE type; - - // inner rows - if (rowIndex > 0 && rowIndex < numRows-1) { - for (int i = 0; i < length; i++) { - sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i)) - * concentrations(rowIndex+1,i) - + ( - 1 - sy * ( - calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,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) - ; - } + // inner rows + if (rowIndex > 0 && rowIndex < numRows - 1) { + for (int i = 0; i < length; i++) { + sv(i) = + sy * + calcAlphaIntercell(alphaY(rowIndex, i), alphaY(rowIndex + 1, i)) * + concentrations(rowIndex + 1, i) + + (1 - sy * (calcAlphaIntercell(alphaY(rowIndex, i), + alphaY(rowIndex + 1, 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); } + } - // first row - else if (rowIndex == 0) { - for (int i = 0; i < length; 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!"); - } - } + // first row + else if (rowIndex == 0) { + for (int i = 0; i < length; 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 - else if (rowIndex == numRows-1) { - for (int i = 0; i < length; 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!"); - } - } + // last row + else if (rowIndex == numRows - 1) { + for (int i = 0; i < length; 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 in constant bc case - if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) { - sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft[rowIndex].getValue(); - } + // 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 in constant bc case - if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) { - sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight[rowIndex].getValue(); - } + // 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; + return sv; } - // solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector // use of EigenLU solver static VectorXd EigenLUAlgorithm(SparseMatrix &A, VectorXd &b) { - SparseLU> solver; - solver.analyzePattern(A); - solver.factorize(A); + SparseLU> solver; + solver.analyzePattern(A); + solver.factorize(A); - return solver.solve(b); + return solver.solve(b); } -// solver for linear equation system; A corresponds to coefficient matrix, +// solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector // implementation of Thomas Algorithm static VectorXd ThomasAlgorithm(SparseMatrix &A, VectorXd &b) { - uint32_t n = b.size(); + 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; + 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); + // 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); + 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]; + // 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]); - } + 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]); + 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]; - } + for (int i = n; i-- > 0;) { + x_vec[i] -= c_diag[i] * x_vec[i + 1]; + } - return x_vec; + return x_vec; } +// BTCS solution for 1D grid +static void BTCS_1D(Grid &grid, Boundary &bc, double timestep, + VectorXd (*solverFunc)(SparseMatrix &A, + VectorXd &b)) { + int length = grid.getLength(); + double sx = timestep / (grid.getDelta() * grid.getDelta()); -// BTCS solution for 1D grid -static void BTCS_1D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solverFunc) (SparseMatrix &A, VectorXd &b)) { - int length = grid.getLength(); - double sx = timestep / (grid.getDelta() * grid.getDelta()); + VectorXd concentrations_t1(length); - VectorXd concentrations_t1(length); + SparseMatrix A; + VectorXd b(length); - SparseMatrix A; - VectorXd b(length); + MatrixXd alpha = grid.getAlpha(); + vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); - MatrixXd alpha = grid.getAlpha(); - vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + MatrixXd concentrations = grid.getConcentrations(); + int rowIndex = 0; + A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, + sx); // this is exactly same as in 2D + for (int i = 0; i < length; i++) { + 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(); + } - MatrixXd concentrations = grid.getConcentrations(); - int rowIndex = 0; - A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx); // this is exactly same as in 2D - for (int i = 0; i < length; i++) { - 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(); - } + concentrations_t1 = solverFunc(A, b); - concentrations_t1 = solverFunc(A, b); + for (int j = 0; j < length; j++) { + concentrations(0, j) = concentrations_t1(j); + } - for (int j = 0; j < length; j++) { - concentrations(0,j) = concentrations_t1(j); - } - - grid.setConcentrations(concentrations); + grid.setConcentrations(concentrations); } - // BTCS solution for 2D grid -static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solverFunc) (SparseMatrix &A, VectorXd &b), int numThreads) { - int rowMax = grid.getRow(); - int colMax = grid.getCol(); - double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); - double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); +static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, + VectorXd (*solverFunc)(SparseMatrix &A, + VectorXd &b), + int numThreads) { + int rowMax = grid.getRow(); + int colMax = grid.getCol(); + double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); + double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); - MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); - VectorXd row_t1(colMax); + MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); + VectorXd row_t1(colMax); - SparseMatrix A; - VectorXd b; + SparseMatrix A; + VectorXd b; - MatrixXd alphaX = grid.getAlphaX(); - MatrixXd alphaY = grid.getAlphaY(); - vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); - vector bcTop = bc.getBoundarySide(BC_SIDE_TOP); - vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); + MatrixXd alphaX = grid.getAlphaX(); + MatrixXd alphaY = grid.getAlphaY(); + vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + vector bcTop = bc.getBoundarySide(BC_SIDE_TOP); + vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); - MatrixXd concentrations = grid.getConcentrations(); + MatrixXd concentrations = grid.getConcentrations(); - #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) - for (int i = 0; i < rowMax; i++) { - - - A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx); - b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, - bcTop, bcBottom, colMax, i, sx, sy); - - SparseLU> solver; +#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) + for (int i = 0; i < rowMax; i++) { - row_t1 = solverFunc(A, b); - - concentrations_t1.row(i) = row_t1; - } - - concentrations_t1.transposeInPlace(); - concentrations.transposeInPlace(); - alphaX.transposeInPlace(); - alphaY.transposeInPlace(); - - #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) - for (int i = 0; i < colMax; i++) { - // swap alphas, boundary conditions and sx/sy for column-wise calculation - A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy); - b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, - bcLeft, bcRight, rowMax, i, sy, sx); - - row_t1 = solverFunc(A, b); + A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx); + b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, + bcTop, bcBottom, colMax, i, sx, sy); - concentrations.row(i) = row_t1; - } - - concentrations.transposeInPlace(); + SparseLU> solver; - grid.setConcentrations(concentrations); + row_t1 = solverFunc(A, b); + + concentrations_t1.row(i) = row_t1; + } + + concentrations_t1.transposeInPlace(); + concentrations.transposeInPlace(); + alphaX.transposeInPlace(); + alphaY.transposeInPlace(); + +#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) + for (int i = 0; i < colMax; i++) { + // swap alphas, boundary conditions and sx/sy for column-wise calculation + A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy); + b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, + bcLeft, bcRight, rowMax, i, sy, sx); + + row_t1 = solverFunc(A, b); + + concentrations.row(i) = row_t1; + } + + concentrations.transposeInPlace(); + + grid.setConcentrations(concentrations); } - // entry point for EigenLU solver; differentiate between 1D and 2D grid static void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { - if (grid.getDim() == 1) { - BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); - } else if (grid.getDim() == 2) { - BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads); - } else { - throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); - } + if (grid.getDim() == 1) { + BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); + } else if (grid.getDim() == 2) { + BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads); + } 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, int numThreads) { - if (grid.getDim() == 1) { - BTCS_1D(grid, bc, timestep, ThomasAlgorithm); - } else if (grid.getDim() == 2) { - BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads); - } else { - throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); - } +static void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, + int numThreads) { + if (grid.getDim() == 1) { + BTCS_1D(grid, bc, timestep, ThomasAlgorithm); + } else if (grid.getDim() == 2) { + BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads); + } else { + throw_invalid_argument( + "Error: Only 1- and 2-dimensional grids are defined!"); + } } \ No newline at end of file diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 2de5ba6..84b81e0 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -1,162 +1,162 @@ #include "TugUtils.cpp" #include #include -#include #include +#include using namespace std; BoundaryElement::BoundaryElement() { - - this->type = BC_TYPE_CLOSED; - this->value = -1; // without meaning in closed case + + this->type = BC_TYPE_CLOSED; + this->value = -1; // without meaning in closed case } BoundaryElement::BoundaryElement(double value) { - this->type = BC_TYPE_CONSTANT; - this->value = value; + this->type = BC_TYPE_CONSTANT; + this->value = value; } -void BoundaryElement::setType(BC_TYPE type) { - this->type = type; -} +void BoundaryElement::setType(BC_TYPE type) { this->type = type; } void BoundaryElement::setValue(double value) { - if(value < 0){ - throw_invalid_argument("No negative concentration allowed."); - } - if(type == BC_TYPE_CLOSED){ - throw_invalid_argument( - "No constant boundary concentrations can be set for closed " - "boundaries. Please change type first."); - } - this->value = value; + if (value < 0) { + throw_invalid_argument("No negative concentration allowed."); + } + if (type == BC_TYPE_CLOSED) { + throw_invalid_argument( + "No constant boundary concentrations can be set for closed " + "boundaries. Please change type first."); + } + this->value = value; } -BC_TYPE BoundaryElement::getType() { - return this->type; -} +BC_TYPE BoundaryElement::getType() { return this->type; } -double BoundaryElement::getValue() { - return this->value; -} +double BoundaryElement::getValue() { return this->value; } Boundary::Boundary(Grid grid) : grid(grid) { - if (grid.getDim() == 1) { - this->boundaries = vector>(2); // in 1D only left and right boundary + if (grid.getDim() == 1) { + this->boundaries = vector>( + 2); // in 1D only left and right boundary - this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); - this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); - } else if (grid.getDim() == 2) { - this->boundaries = vector>(4); + this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); + } else if (grid.getDim() == 2) { + this->boundaries = vector>(4); - this->boundaries[BC_SIDE_LEFT] = vector(grid.getRow(), BoundaryElement()); - this->boundaries[BC_SIDE_RIGHT] = vector(grid.getRow(), BoundaryElement()); - this->boundaries[BC_SIDE_TOP] = vector(grid.getCol(), BoundaryElement()); - this->boundaries[BC_SIDE_BOTTOM] = vector(grid.getCol(), BoundaryElement()); - } + this->boundaries[BC_SIDE_LEFT] = + vector(grid.getRow(), BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT] = + vector(grid.getRow(), BoundaryElement()); + this->boundaries[BC_SIDE_TOP] = + vector(grid.getCol(), BoundaryElement()); + this->boundaries[BC_SIDE_BOTTOM] = + vector(grid.getCol(), BoundaryElement()); + } } void Boundary::setBoundarySideClosed(BC_SIDE side) { - if(grid.getDim() == 1){ - if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ - throw_invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } + if (grid.getDim() == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw_invalid_argument( + "For the one-dimensional case, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); } + } - int n; - if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) { - n = grid.getRow(); - } else { - n = grid.getCol(); - } - this->boundaries[side] = vector(n, BoundaryElement()); + int n; + if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) { + n = grid.getRow(); + } else { + n = grid.getCol(); + } + this->boundaries[side] = vector(n, BoundaryElement()); } void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { - if(grid.getDim() == 1){ - if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ - throw_invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } + if (grid.getDim() == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw_invalid_argument( + "For the one-dimensional case, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); } + } - int n; - if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) { - n = grid.getRow(); - } else { - n = grid.getCol(); - } - this->boundaries[side] = vector(n, BoundaryElement(value)); + int n; + if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) { + n = grid.getRow(); + } else { + n = grid.getCol(); + } + this->boundaries[side] = vector(n, BoundaryElement(value)); } void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) { - // tests whether the index really points to an element of the boundary side. - if((boundaries[side].size() < index) || index < 0){ - throw_invalid_argument("Index is selected either too large or too small."); - } - this->boundaries[side][index].setType(BC_TYPE_CLOSED); + // tests whether the index really points to an element of the boundary side. + if ((boundaries[side].size() < index) || index < 0) { + throw_invalid_argument("Index is selected either too large or too small."); + } + this->boundaries[side][index].setType(BC_TYPE_CLOSED); } -void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, double value) { - // tests whether the index really points to an element of the boundary side. - if((boundaries[side].size() < index) || index < 0){ - throw_invalid_argument("Index is selected either too large or too small."); - } - this->boundaries[side][index].setType(BC_TYPE_CONSTANT); - this->boundaries[side][index].setValue(value); +void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, + double value) { + // tests whether the index really points to an element of the boundary side. + if ((boundaries[side].size() < index) || index < 0) { + throw_invalid_argument("Index is selected either too large or too small."); + } + this->boundaries[side][index].setType(BC_TYPE_CONSTANT); + this->boundaries[side][index].setValue(value); } const vector Boundary::getBoundarySide(BC_SIDE side) { - if(grid.getDim() == 1){ - if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ - throw_invalid_argument( - "For the one-dimensional trap, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } + if (grid.getDim() == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw_invalid_argument( + "For the one-dimensional trap, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); } - return this->boundaries[side]; + } + return this->boundaries[side]; } VectorXd Boundary::getBoundarySideValues(BC_SIDE side) { - int length = boundaries[side].size(); - VectorXd values(length); + int length = boundaries[side].size(); + VectorXd values(length); - for (int i = 0; i < length; i++) { - if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { - values(i) = -1; - continue; - } - values(i) = getBoundaryElementValue(side, i); + for (int i = 0; i < length; i++) { + if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { + values(i) = -1; + continue; } + values(i) = getBoundaryElementValue(side, i); + } - return values; + return values; } BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) { - if((boundaries[side].size() < index) || index < 0){ - throw_invalid_argument("Index is selected either too large or too small."); - } - return this->boundaries[side][index]; + if ((boundaries[side].size() < index) || index < 0) { + throw_invalid_argument("Index is selected either too large or too small."); + } + return this->boundaries[side][index]; } BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) { - if((boundaries[side].size() < index) || index < 0){ - throw_invalid_argument("Index is selected either too large or too small."); - } - return this->boundaries[side][index].getType(); + if ((boundaries[side].size() < index) || index < 0) { + throw_invalid_argument("Index is selected either too large or too small."); + } + return this->boundaries[side][index].getType(); } double Boundary::getBoundaryElementValue(BC_SIDE side, int index) { - if((boundaries[side].size() < index) || index < 0){ - throw_invalid_argument("Index is selected either too large or too small."); - } - if(boundaries[side][index].getType() != BC_TYPE_CONSTANT){ - throw_invalid_argument("A value can only be output if it is a constant boundary condition."); - } - return this->boundaries[side][index].getValue(); + if ((boundaries[side].size() < index) || index < 0) { + throw_invalid_argument("Index is selected either too large or too small."); + } + if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { + throw_invalid_argument( + "A value can only be output if it is a constant boundary condition."); + } + return this->boundaries[side][index].getValue(); } - diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 647c885..d48cc46 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -1,433 +1,389 @@ /** * @file FTCS.cpp - * @brief Implementation of heterogenous FTCS (forward time-centered space) solution - * of diffusion equation in 1D and 2D space. - * + * @brief Implementation of heterogenous FTCS (forward time-centered space) + * solution of diffusion equation in 1D and 2D space. + * */ #include "TugUtils.cpp" #include -#include #include #include +#include using namespace std; - // calculates horizontal change on one cell independent of boundary type static double calcHorizontalChange(Grid &grid, int &row, int &col) { - double result = - calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col+1) - - ( - calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) - ) - * grid.getConcentrations()(row,col) - + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col-1); + double result = calcAlphaIntercell(grid.getAlphaX()(row, col + 1), + grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col + 1) - + (calcAlphaIntercell(grid.getAlphaX()(row, col + 1), + grid.getAlphaX()(row, col)) + + calcAlphaIntercell(grid.getAlphaX()(row, col - 1), + grid.getAlphaX()(row, col))) * + grid.getConcentrations()(row, col) + + calcAlphaIntercell(grid.getAlphaX()(row, col - 1), + grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col - 1); - return result; + return result; } - // calculates vertical change on one cell independent of boundary type static double calcVerticalChange(Grid &grid, int &row, int &col) { - - double result = - calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) - * grid.getConcentrations()(row+1,col) - - ( - calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) - + calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)) - ) - * grid.getConcentrations()(row,col) - + calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)) - * grid.getConcentrations()(row-1,col); - return result; + double result = calcAlphaIntercell(grid.getAlphaY()(row + 1, col), + grid.getAlphaY()(row, col)) * + grid.getConcentrations()(row + 1, col) - + (calcAlphaIntercell(grid.getAlphaY()(row + 1, col), + grid.getAlphaY()(row, col)) + + calcAlphaIntercell(grid.getAlphaY()(row - 1, col), + grid.getAlphaY()(row, col))) * + grid.getConcentrations()(row, col) + + calcAlphaIntercell(grid.getAlphaY()(row - 1, col), + grid.getAlphaY()(row, col)) * + grid.getConcentrations()(row - 1, col); + + return result; } - // calculates horizontal change on one cell with a constant left boundary -static double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { +static double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, Boundary &bc, + int &row, int &col) { - double result = - calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col+1) - - ( - calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) - + 2 * grid.getAlphaX()(row,col) - ) - * grid.getConcentrations()(row,col) - + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryElementValue(BC_SIDE_LEFT, row); + double result = calcAlphaIntercell(grid.getAlphaX()(row, col + 1), + grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col + 1) - + (calcAlphaIntercell(grid.getAlphaX()(row, col + 1), + grid.getAlphaX()(row, col)) + + 2 * grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col) + + 2 * grid.getAlphaX()(row, col) * + bc.getBoundaryElementValue(BC_SIDE_LEFT, row); - return result; + return result; } - // calculates horizontal change on one cell with a closed left boundary -static 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; -} +static 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; +} // checks boundary condition type for a cell on the left edge of grid -static 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!"); - } +static 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!"); + } } - // calculates horizontal change on one cell with a constant right boundary -static double calcHorizontalChangeRightBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { +static double calcHorizontalChangeRightBoundaryConstant(Grid &grid, + Boundary &bc, int &row, + int &col) { - double result = - 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) - ) - * grid.getConcentrations()(row,col) - + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)) - * grid.getConcentrations()(row,col-1); + double result = 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)) * + grid.getConcentrations()(row, col) + + calcAlphaIntercell(grid.getAlphaX()(row, col - 1), + grid.getAlphaX()(row, col)) * + grid.getConcentrations()(row, col - 1); - return result; + return result; } - // calculates horizontal change on one cell with a closed right boundary -static 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; -} +static 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; +} // checks boundary condition type for a cell on the right edge of grid -static 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!"); - } +static 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!"); + } } - // calculates vertical change on one cell with a constant top boundary -static double calcVerticalChangeTopBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { - - double result = - calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col)) - * grid.getConcentrations()(row+1,col) - - ( - calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col)) - + 2 * grid.getAlphaY()(row, col) - ) - * grid.getConcentrations()(row, col) - + 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_TOP, col); +static double calcVerticalChangeTopBoundaryConstant(Grid &grid, Boundary &bc, + int &row, int &col) { - return result; + double result = calcAlphaIntercell(grid.getAlphaY()(row + 1, col), + grid.getAlphaY()(row, col)) * + grid.getConcentrations()(row + 1, col) - + (calcAlphaIntercell(grid.getAlphaY()(row + 1, col), + grid.getAlphaY()(row, col)) + + 2 * grid.getAlphaY()(row, col)) * + grid.getConcentrations()(row, col) + + 2 * grid.getAlphaY()(row, col) * + bc.getBoundaryElementValue(BC_SIDE_TOP, col); + + return result; } - // calculates vertical change on one cell with a closed top boundary -static 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)); +static double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, + int &col) { - return result; + double result = calcAlphaIntercell(grid.getAlphaY()(row + 1, col), + grid.getConcentrations()(row, col)) * + (grid.getConcentrations()(row + 1, col) - + grid.getConcentrations()(row, col)); + + return result; } - // checks boundary condition type for a cell on the top edge of grid -static 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!"); - } +static 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!"); + } } - // calculates vertical change on one cell with a constant bottom boundary -static double calcVerticalChangeBottomBoundaryConstant(Grid &grid, Boundary &bc, int &row, int &col) { +static double calcVerticalChangeBottomBoundaryConstant(Grid &grid, Boundary &bc, + int &row, int &col) { - double result = - 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) - ) - * grid.getConcentrations()(row, col) - + calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) - * grid.getConcentrations()(row-1,col); + double result = 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)) * + grid.getConcentrations()(row, col) + + calcAlphaIntercell(grid.getAlphaY()(row, col), + grid.getAlphaY()(row - 1, col)) * + grid.getConcentrations()(row - 1, col); - return result; + return result; } - // calculates vertical change on one cell with a closed bottom boundary -static double calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) { +static 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))); + double result = -(calcAlphaIntercell(grid.getAlphaY()(row, col), + grid.getAlphaY()(row - 1, col)) * + (grid.getConcentrations()(row, col) - + grid.getConcentrations()(row - 1, col))); - return result; + return result; } - // checks boundary condition type for a cell on the bottom edge of grid -static 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!"); - } +static 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!"); + } } - // FTCS solution for 1D grid static void FTCS_1D(Grid &grid, Boundary &bc, double ×tep) { - int colMax = grid.getCol(); - double deltaCol = grid.getDeltaCol(); + int colMax = grid.getCol(); + double deltaCol = grid.getDeltaCol(); - // matrix for concentrations at time t+1 - MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0); + // matrix for concentrations at time t+1 + MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0); - // only one row in 1D case -> row constant at index 0 - int row = 0; + // only one row in 1D case -> row constant at index 0 + int row = 0; - // inner cells - // independent of boundary condition type - for (int col = 1; col < colMax-1; col++) { - concentrations_t1(row,col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChange(grid, row, col) - ) - ; - } + // inner cells + // independent of boundary condition type + for (int col = 1; col < colMax - 1; col++) { + concentrations_t1(row, col) = grid.getConcentrations()(row, col) + + timestep / (deltaCol * deltaCol) * + (calcHorizontalChange(grid, row, col)); + } - // left boundary; hold column constant at index 0 - int col = 0; - concentrations_t1(row, col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChangeLeftBoundary(grid, bc, row, col) - ) - ; + // left boundary; hold column constant at index 0 + int col = 0; + concentrations_t1(row, col) = + grid.getConcentrations()(row, col) + + timestep / (deltaCol * deltaCol) * + (calcHorizontalChangeLeftBoundary(grid, bc, row, col)); + // right boundary; hold column constant at max index + col = colMax - 1; + concentrations_t1(row, col) = + grid.getConcentrations()(row, col) + + timestep / (deltaCol * deltaCol) * + (calcHorizontalChangeRightBoundary(grid, bc, row, col)); - // right boundary; hold column constant at max index - col = colMax-1; - concentrations_t1(row,col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChangeRightBoundary(grid, bc, row, col) - ) - ; - - // overwrite obsolete concentrations - grid.setConcentrations(concentrations_t1); + // overwrite obsolete concentrations + grid.setConcentrations(concentrations_t1); } - // FTCS solution for 2D grid -static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep, int numThreads) { - int rowMax = grid.getRow(); - int colMax = grid.getCol(); - double deltaRow = grid.getDeltaRow(); - double deltaCol = grid.getDeltaCol(); +static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep, + int numThreads) { + int rowMax = grid.getRow(); + int colMax = grid.getCol(); + double deltaRow = grid.getDeltaRow(); + double deltaCol = grid.getDeltaCol(); - // matrix for concentrations at time t+1 - MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); - - // inner cells - // these are independent of the boundary condition type - // omp_set_num_threads(10); -#pragma omp parallel for num_threads(numThreads) - for (int row = 1; row < rowMax-1; row++) { - for (int col = 1; col < colMax-1; col++) { - concentrations_t1(row, col) = grid.getConcentrations()(row, col) - + timestep / (deltaRow*deltaRow) - * ( - calcVerticalChange(grid, row, col) - ) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChange(grid, row, col) - ) - ; - } - } - - // boundary conditions - // left without corners / looping over rows - // hold column constant at index 0 - int col = 0; -#pragma omp parallel for num_threads(numThreads) - for (int row = 1; row < rowMax-1; row++) { - concentrations_t1(row, col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChangeLeftBoundary(grid, bc, row, col) - ) - + timestep / (deltaRow*deltaRow) - * ( - calcVerticalChange(grid, row, col) - ) - ; - } - - // right without corners / looping over rows - // hold column constant at max index - col = colMax-1; -#pragma omp parallel for num_threads(numThreads) - for (int row = 1; row < rowMax-1; row++) { - concentrations_t1(row,col) = grid.getConcentrations()(row,col) - + timestep / (deltaCol*deltaCol) - * ( - calcHorizontalChangeRightBoundary(grid, bc, row, col) - ) - + timestep / (deltaRow*deltaRow) - * ( - calcVerticalChange(grid, row, col) - ) - ; - } - - - // top without corners / looping over columns - // hold row constant at index 0 - int row = 0; -#pragma omp parallel for num_threads(numThreads) - for (int col=1; col #include +#include Grid::Grid(int length) { - if (length <= 3) { - throw_invalid_argument("Given grid length too small. Must be greater than 3."); - } + if (length <= 3) { + throw_invalid_argument( + "Given grid length too small. Must be greater than 3."); + } - this->row = 1; - this->col = length; - this->domainCol = length; // default: same size as length - this->deltaCol = double(this->domainCol)/double(this->col); // -> 1 - this->dim = 1; + this->row = 1; + this->col = length; + this->domainCol = length; // default: same size as length + this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + this->dim = 1; - this->concentrations = MatrixXd::Constant(1, col, 20); - this->alphaX = MatrixXd::Constant(1, col, 1); + this->concentrations = MatrixXd::Constant(1, col, 20); + this->alphaX = MatrixXd::Constant(1, col, 1); } Grid::Grid(int row, int col) { - if (row <= 3 || col <= 3) { - throw_invalid_argument("Given grid dimensions too small. Must each be greater than 3."); - } + if (row <= 3 || col <= 3) { + throw_invalid_argument( + "Given grid dimensions too small. Must each be greater than 3."); + } - this->row = row; - this->col = col; - this->domainRow = row; // default: same size as row - this->domainCol = col; // default: same size as col - this->deltaRow = double(this->domainRow)/double(this->row); // -> 1 - this->deltaCol = double(this->domainCol)/double(this->col); // -> 1 - this->dim = 2; + this->row = row; + this->col = col; + this->domainRow = row; // default: same size as row + this->domainCol = col; // default: same size as col + this->deltaRow = double(this->domainRow) / double(this->row); // -> 1 + this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + this->dim = 2; - this->concentrations = MatrixXd::Constant(row, col, 20); - this->alphaX = MatrixXd::Constant(row, col, 1); - this->alphaY = MatrixXd::Constant(row, col, 1); + this->concentrations = MatrixXd::Constant(row, col, 20); + this->alphaX = MatrixXd::Constant(row, col, 1); + this->alphaY = MatrixXd::Constant(row, col, 1); } void Grid::setConcentrations(MatrixXd concentrations) { - if (concentrations.rows() != this->row || concentrations.cols() != this->col) { - throw_invalid_argument("Given matrix of concentrations mismatch with Grid dimensions!"); - } + if (concentrations.rows() != this->row || + concentrations.cols() != this->col) { + throw_invalid_argument( + "Given matrix of concentrations mismatch with Grid dimensions!"); + } - this->concentrations = concentrations; + this->concentrations = concentrations; } -const MatrixXd Grid::getConcentrations() { - return this->concentrations; -} +const MatrixXd Grid::getConcentrations() { return this->concentrations; } void Grid::setAlpha(MatrixXd alpha) { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably use 2D setter function!"); - } - if (alpha.rows() != 1 || alpha.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients mismatch with Grid dimensions!"); - } + if (dim != 1) { + throw_invalid_argument("Grid is not one dimensional, you should probably " + "use 2D setter function!"); + } + if (alpha.rows() != 1 || alpha.cols() != this->col) { + throw_invalid_argument( + "Given matrix of alpha coefficients mismatch with Grid dimensions!"); + } - this->alphaX = alpha; + this->alphaX = alpha; } void Grid::setAlpha(MatrixXd alphaX, MatrixXd alphaY) { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably use 1D setter function!"); - } - if (alphaX.rows() != this->row || alphaX.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients in x-direction mismatch with GRid dimensions!"); - } - if (alphaY.rows() != this->row || alphaY.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients in y-direction mismatch with GRid dimensions!"); - } + if (dim != 2) { + throw_invalid_argument("Grid is not two dimensional, you should probably " + "use 1D setter function!"); + } + if (alphaX.rows() != this->row || alphaX.cols() != this->col) { + throw_invalid_argument("Given matrix of alpha coefficients in x-direction " + "mismatch with GRid dimensions!"); + } + if (alphaY.rows() != this->row || alphaY.cols() != this->col) { + throw_invalid_argument("Given matrix of alpha coefficients in y-direction " + "mismatch with GRid dimensions!"); + } - this->alphaX = alphaX; - this->alphaY = alphaY; + this->alphaX = alphaX; + this->alphaY = alphaY; } const MatrixXd Grid::getAlpha() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably use either getAlphaX() or getAlphaY()!"); - } + if (dim != 1) { + throw_invalid_argument("Grid is not one dimensional, you should probably " + "use either getAlphaX() or getAlphaY()!"); + } - return this->alphaX; + return this->alphaX; } const MatrixXd Grid::getAlphaX() { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably use getAlpha()!"); - } + if (dim != 2) { + throw_invalid_argument( + "Grid is not two dimensional, you should probably use getAlpha()!"); + } - return this->alphaX; + return this->alphaX; } const MatrixXd Grid::getAlphaY() { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably use getAlpha()!"); - } + if (dim != 2) { + throw_invalid_argument( + "Grid is not two dimensional, you should probably use getAlpha()!"); + } - return this->alphaY; + return this->alphaY; } -int Grid::getDim() { - return dim; -} +int Grid::getDim() { return dim; } int Grid::getLength() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably use getRow() or getCol()!"); - } + if (dim != 1) { + throw_invalid_argument("Grid is not one dimensional, you should probably " + "use getRow() or getCol()!"); + } - return col; + return col; } -int Grid::getRow() { - return row; -} +int Grid::getRow() { return row; } -int Grid::getCol() { - return col; -} +int Grid::getCol() { return col; } void Grid::setDomain(double domainLength) { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probaly use the 2D domain setter!"); - } - if (domainLength <= 0) { - throw_invalid_argument("Given domain length is not positive!"); - } + if (dim != 1) { + throw_invalid_argument("Grid is not one dimensional, you should probaly " + "use the 2D domain setter!"); + } + if (domainLength <= 0) { + throw_invalid_argument("Given domain length is not positive!"); + } - this->domainCol = domainLength; - this->deltaCol = double(this->domainCol)/double(this->col); + this->domainCol = domainLength; + this->deltaCol = double(this->domainCol) / double(this->col); } void Grid::setDomain(double domainRow, double domainCol) { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably use the 1D domain setter!"); - } - if (domainRow <= 0 || domainCol <= 0) { - throw_invalid_argument("Given domain size is not positive!"); - } + if (dim != 2) { + throw_invalid_argument("Grid is not two dimensional, you should probably " + "use the 1D domain setter!"); + } + if (domainRow <= 0 || domainCol <= 0) { + throw_invalid_argument("Given domain size is not positive!"); + } - this->domainRow = domainRow; - this->domainCol = domainCol; - this->deltaRow = double(this->domainRow)/double(this->row); - this->deltaCol = double(this->domainCol)/double(this->col); + this->domainRow = domainRow; + this->domainCol = domainCol; + this->deltaRow = double(this->domainRow) / double(this->row); + 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"); - } + if (dim != 1) { + throw_invalid_argument("Grid is not one dimensional, you should probably " + "use the 2D delta getters"); + } - return this->deltaCol; + return this->deltaCol; } -double Grid::getDeltaCol() { - return this->deltaCol; -} +double Grid::getDeltaCol() { return this->deltaCol; } double Grid::getDeltaRow() { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, meaning there is no delta in y-direction!"); - } + if (dim != 2) { + throw_invalid_argument("Grid is not two dimensional, meaning there is no " + "delta in y-direction!"); + } - return this->deltaRow; + return this->deltaRow; } diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 09e5198..31d1724 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -1,10 +1,10 @@ #include #include #include +#include #include #include #include -#include #include @@ -13,319 +13,330 @@ #include "BTCSv2.cpp" - using namespace std; -Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid), bc(bc) { +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; needs to be set - this->iterations = -1; - this->innerIterations = 1; - this->numThreads = omp_get_num_procs()-1; - - this->csv_output = CSV_OUTPUT_OFF; - this->console_output = CONSOLE_OUTPUT_OFF; - this->time_measure = TIME_MEASURE_OFF; + this->approach = approach; + this->solver = THOMAS_ALGORITHM_SOLVER; + this->timestep = -1; // error per default; needs to be set + this->iterations = -1; + this->innerIterations = 1; + this->numThreads = omp_get_num_procs() - 1; + + this->csv_output = CSV_OUTPUT_OFF; + this->console_output = CONSOLE_OUTPUT_OFF; + this->time_measure = TIME_MEASURE_OFF; } 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!"); - } + if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { + throw_invalid_argument("Invalid CSV output option given!"); + } - this->csv_output = csv_output; + this->csv_output = csv_output; } void Simulation::setOutputConsole(CONSOLE_OUTPUT console_output) { - if (console_output < CONSOLE_OUTPUT_OFF && console_output > CONSOLE_OUTPUT_VERBOSE) { - throw_invalid_argument("Invalid console output option given!"); - } + if (console_output < CONSOLE_OUTPUT_OFF && + console_output > CONSOLE_OUTPUT_VERBOSE) { + throw_invalid_argument("Invalid console output option given!"); + } - this->console_output = console_output; + this->console_output = console_output; } void Simulation::setTimeMeasure(TIME_MEASURE time_measure) { - if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { - throw_invalid_argument("Invalid time measure option given!"); - } + if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { + throw_invalid_argument("Invalid time measure option given!"); + } - this->time_measure = time_measure; + this->time_measure = time_measure; } void Simulation::setTimestep(double timestep) { - if(timestep <= 0){ - throw_invalid_argument("Timestep has to be greater than zero."); - } + if (timestep <= 0) { + throw_invalid_argument("Timestep has to be greater than zero."); + } - if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { + if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { - double deltaRowSquare; - double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - double minDeltaSquare; - double maxAlphaX, maxAlphaY, maxAlpha; - string dim; - if (grid.getDim() == 2) { - dim = "2D"; + double deltaRowSquare; + double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + double minDeltaSquare; + double maxAlphaX, maxAlphaY, maxAlpha; + string dim; + if (grid.getDim() == 2) { + dim = "2D"; - deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); - minDeltaSquare = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - maxAlphaX = grid.getAlphaX().maxCoeff(); - maxAlphaY = grid.getAlphaY().maxCoeff(); - maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; + minDeltaSquare = + (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; + maxAlphaX = grid.getAlphaX().maxCoeff(); + maxAlphaY = grid.getAlphaY().maxCoeff(); + maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; - } else if (grid.getDim() == 1) { - dim = "1D"; - minDeltaSquare = deltaColSquare; - maxAlpha = grid.getAlpha().maxCoeff(); - - } else { - throw_invalid_argument("Critical error: Undefined number of dimensions!"); - } - - // Courant-Friedrichs-Lewy condition - double cfl = minDeltaSquare / (4*maxAlpha); - - // 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))); - - string approachPrefix = (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI"); - cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl << endl; - cout << approachPrefix << "_" << dim << " :: required dt=" << timestep << endl; - - if (timestep > cfl) { - - 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 << approachPrefix << "_" << dim << " :: Required " << this->innerIterations - << " inner iterations with dt=" << this->timestep << endl; - - } else { - - this->timestep = timestep; - cout << approachPrefix << "_" << dim << " :: No inner iterations required, dt=" << timestep << endl; - - } + } else if (grid.getDim() == 1) { + dim = "1D"; + minDeltaSquare = deltaColSquare; + maxAlpha = grid.getAlpha().maxCoeff(); } else { - this->timestep = timestep; + throw_invalid_argument("Critical error: Undefined number of dimensions!"); } + // Courant-Friedrichs-Lewy condition + double cfl = minDeltaSquare / (4 * maxAlpha); + + // 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))); + + string approachPrefix = + (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI"); + cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl + << endl; + cout << approachPrefix << "_" << dim << " :: required dt=" << timestep + << endl; + + if (timestep > cfl) { + + 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 << approachPrefix << "_" << dim << " :: Required " + << this->innerIterations + << " inner iterations with dt=" << this->timestep << endl; + + } else { + + this->timestep = timestep; + cout << approachPrefix << "_" << dim + << " :: No inner iterations required, dt=" << timestep << endl; + } + + } else { + this->timestep = timestep; + } } -double Simulation::getTimestep() { - return this->timestep; -} +double Simulation::getTimestep() { return this->timestep; } void Simulation::setIterations(int iterations) { - if(iterations <= 0){ - throw_invalid_argument("Number of iterations must be greater than zero."); - } - this->iterations = iterations; + if (iterations <= 0) { + throw_invalid_argument("Number of iterations must be greater than zero."); + } + 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; - } + 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; + this->solver = solver; } -void Simulation::setNumberThreads(int numThreads){ - if(numThreads>0 && numThreads<=omp_get_num_procs()){ - this->numThreads=numThreads; - } - else{ - int maxThreadNumber = omp_get_num_procs(); - string outputMessage = "Number of threads exceeds the number of processor cores (" - + to_string(maxThreadNumber) + ") or is less than 1."; - - throw_invalid_argument(outputMessage); - } +void Simulation::setNumberThreads(int numThreads) { + if (numThreads > 0 && numThreads <= omp_get_num_procs()) { + this->numThreads = numThreads; + } else { + int maxThreadNumber = omp_get_num_procs(); + string outputMessage = + "Number of threads exceeds the number of processor cores (" + + to_string(maxThreadNumber) + ") or is less than 1."; + + throw_invalid_argument(outputMessage); + } } -int Simulation::getIterations() { - return this->iterations; -} +int Simulation::getIterations() { return this->iterations; } void Simulation::printConcentrationsConsole() { - cout << grid.getConcentrations() << endl; - cout << endl; + cout << grid.getConcentrations() << endl; + cout << endl; } string Simulation::createCSVfile() { - ofstream file; - int appendIdent = 0; - string appendIdentString; + ofstream file; + int appendIdent = 0; + string appendIdentString; - // string approachString = (approach == 0) ? "FTCS" : "BTCS"; - string approachString = (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI"); - string row = to_string(grid.getRow()); - string col = to_string(grid.getCol()); - string numIterations = to_string(iterations); + // string approachString = (approach == 0) ? "FTCS" : "BTCS"; + string approachString = + (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI"); + string row = to_string(grid.getRow()); + string col = to_string(grid.getCol()); + string numIterations = to_string(iterations); - string filename = approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; + string filename = + approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; - while (filesystem::exists(filename)) { - appendIdent += 1; - appendIdentString = to_string(appendIdent); - filename = approachString + "_" + row + "_" + col + "_" + numIterations + "-" + appendIdentString + ".csv"; - } + while (filesystem::exists(filename)) { + appendIdent += 1; + appendIdentString = to_string(appendIdent); + filename = approachString + "_" + row + "_" + col + "_" + numIterations + + "-" + appendIdentString + ".csv"; + } - file.open(filename); - if (!file) { - exit(1); - } + file.open(filename); + if (!file) { + 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) { - 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; - } + // 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) { + 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; + } - file.close(); + file.close(); - return filename; + return filename; } void Simulation::printConcentrationsCSV(string filename) { - ofstream file; + ofstream file; - file.open(filename, std::ios_base::app); - if (!file) { - exit(1); - } + file.open(filename, std::ios_base::app); + if (!file) { + exit(1); + } - IOFormat do_not_align(StreamPrecision, DontAlignCols); - file << grid.getConcentrations().format(do_not_align) << endl; - file << endl << endl; - file.close(); + IOFormat do_not_align(StreamPrecision, DontAlignCols); + file << grid.getConcentrations().format(do_not_align) << endl; + file << endl << endl; + file.close(); } void Simulation::run() { - if (this->timestep == -1) { - throw_invalid_argument("Timestep is not set!"); - } - if (this->iterations == -1) { - throw_invalid_argument("Number of iterations are not set!"); - } + if (this->timestep == -1) { + throw_invalid_argument("Timestep is not set!"); + } + if (this->iterations == -1) { + throw_invalid_argument("Number of iterations are not set!"); + } - string filename; - if (this->console_output > CONSOLE_OUTPUT_OFF) { + string filename; + if (this->console_output > CONSOLE_OUTPUT_OFF) { + printConcentrationsConsole(); + } + if (this->csv_output > CSV_OUTPUT_OFF) { + filename = createCSVfile(); + } + + auto begin = std::chrono::high_resolution_clock::now(); + + if (approach == FTCS_APPROACH) { // FTCS case + + for (int i = 0; i < iterations * innerIterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { printConcentrationsConsole(); - } - if (this->csv_output > CSV_OUTPUT_OFF) { - filename = createCSVfile(); - } - - auto begin = std::chrono::high_resolution_clock::now(); - - if (approach == FTCS_APPROACH) { // FTCS 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); - } - - FTCS(this->grid, this->bc, this->timestep, this->numThreads); - - // if (i % (iterations * innerIterations / 100) == 0) { - // double percentage = (double)i / ((double)iterations * (double)innerIterations) * 100; - // if ((int)percentage % 10 == 0) { - // cout << "Progress: " << percentage << "%" << endl; - // } - // } - } - - } else if (approach == BTCS_APPROACH) { // BTCS case - - 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_LU(this->grid, this->bc, this->timestep, this->numThreads); - - } - } 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, this->numThreads); - - } - } - - } else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case - - double beta = 0.5; - - // TODO this implementation is very inefficient! - // a separate implementation that sets up a specific tridiagonal matrix for Crank-Nicolson would be better - MatrixXd concentrations; - MatrixXd concentrationsFTCS; - MatrixXd concentrationsResult; - 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); - } - - concentrations = grid.getConcentrations(); - FTCS(this->grid, this->bc, this->timestep, this->numThreads); - concentrationsFTCS = grid.getConcentrations(); - grid.setConcentrations(concentrations); - BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); - concentrationsResult = beta * concentrationsFTCS + (1-beta) * grid.getConcentrations(); - grid.setConcentrations(concentrationsResult); - } - - } - - auto end = std::chrono::high_resolution_clock::now(); - auto milliseconds = std::chrono::duration_cast(end - begin); - - if (this->console_output > CONSOLE_OUTPUT_OFF) { - printConcentrationsConsole(); - } - if (this->csv_output > CSV_OUTPUT_OFF) { + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { printConcentrationsCSV(filename); + } + + FTCS(this->grid, this->bc, this->timestep, this->numThreads); + + // if (i % (iterations * innerIterations / 100) == 0) { + // double percentage = (double)i / ((double)iterations * + // (double)innerIterations) * 100; if ((int)percentage % 10 == 0) { + // cout << "Progress: " << percentage << "%" << endl; + // } + // } } - if (this->time_measure > TIME_MEASURE_OFF) { - string approachString = (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI"); - string dimString = (grid.getDim() == 1) ? "-1D" : "-2D"; - cout << approachString << dimString << ":: run() finished in " << milliseconds.count() << "ms" << endl; + + } else if (approach == BTCS_APPROACH) { // BTCS case + + 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_LU(this->grid, this->bc, this->timestep, this->numThreads); + } + } 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, this->numThreads); + } } - + + } else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case + + double beta = 0.5; + + // TODO this implementation is very inefficient! + // a separate implementation that sets up a specific tridiagonal matrix for + // Crank-Nicolson would be better + MatrixXd concentrations; + MatrixXd concentrationsFTCS; + MatrixXd concentrationsResult; + 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); + } + + concentrations = grid.getConcentrations(); + FTCS(this->grid, this->bc, this->timestep, this->numThreads); + concentrationsFTCS = grid.getConcentrations(); + grid.setConcentrations(concentrations); + BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); + concentrationsResult = + beta * concentrationsFTCS + (1 - beta) * grid.getConcentrations(); + grid.setConcentrations(concentrationsResult); + } + } + + auto end = std::chrono::high_resolution_clock::now(); + auto milliseconds = + std::chrono::duration_cast(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" : ((approach == 1) ? "BTCS" : "CRNI"); + string dimString = (grid.getDim() == 1) ? "-1D" : "-2D"; + cout << approachString << dimString << ":: run() finished in " + << milliseconds.count() << "ms" << endl; + } } #endif diff --git a/src/TugUtils.cpp b/src/TugUtils.cpp index c75d897..ac0a41e 100644 --- a/src/TugUtils.cpp +++ b/src/TugUtils.cpp @@ -1,7 +1,7 @@ #include +#include #include #include -#include using namespace std; @@ -29,10 +29,11 @@ using namespace std; }) // calculates arithmetic or harmonic mean of alpha between two cells -static double calcAlphaIntercell(const double &alpha1, const double &alpha2, bool useHarmonic = true) { - if (useHarmonic) { - return double(2) / ((double(1)/alpha1) + (double(1)/alpha2)); - } else { - return 0.5 * (alpha1 + alpha2); - } +static double calcAlphaIntercell(const double &alpha1, const double &alpha2, + bool useHarmonic = true) { + if (useHarmonic) { + return double(2) / ((double(1) / alpha1) + (double(1) / alpha2)); + } else { + return 0.5 * (alpha1 + alpha2); + } } diff --git a/test/TestUtils.cpp b/test/TestUtils.cpp index ef0d3db..c6a3662 100644 --- a/test/TestUtils.cpp +++ b/test/TestUtils.cpp @@ -1,8 +1,8 @@ -#include -#include #include #include #include +#include +#include #include #include @@ -11,37 +11,39 @@ using namespace Eigen; MatrixXd CSV2Eigen(string file2Convert) { - vector matrixEntries; + vector matrixEntries; - ifstream matrixDataFile(file2Convert); - if (matrixDataFile.fail()) { - throw invalid_argument("File probably non-existent!"); + ifstream matrixDataFile(file2Convert); + if (matrixDataFile.fail()) { + throw invalid_argument("File probably non-existent!"); + } + + string matrixRowString; + string matrixEntry; + int matrixRowNumber = 0; + + while (getline(matrixDataFile, matrixRowString)) { + stringstream matrixRowStringStream(matrixRowString); + while (getline(matrixRowStringStream, matrixEntry, ' ')) { + matrixEntries.push_back(stod(matrixEntry)); } - - string matrixRowString; - string matrixEntry; - int matrixRowNumber = 0; - - while(getline(matrixDataFile, matrixRowString)){ - stringstream matrixRowStringStream(matrixRowString); - while(getline(matrixRowStringStream, matrixEntry, ' ')){ - matrixEntries.push_back(stod(matrixEntry)); - } - if (matrixRowString.length() > 1) { - matrixRowNumber++; - } + if (matrixRowString.length() > 1) { + matrixRowNumber++; } + } - return Map>(matrixEntries.data(), matrixRowNumber, matrixEntries.size() / matrixRowNumber); + return Map>( + matrixEntries.data(), matrixRowNumber, + matrixEntries.size() / matrixRowNumber); } -bool checkSimilarity(MatrixXd a, MatrixXd b, double precision=1e-5) { - return a.isApprox(b, precision); +bool checkSimilarity(MatrixXd a, MatrixXd b, double precision = 1e-5) { + return a.isApprox(b, precision); } bool checkSimilarityV2(MatrixXd a, MatrixXd b, double maxDiff) { - MatrixXd diff = a - b; - double maxCoeff = diff.maxCoeff(); - return abs(maxCoeff) < maxDiff; + MatrixXd diff = a - b; + double maxCoeff = diff.maxCoeff(); + return abs(maxCoeff) < maxDiff; } \ No newline at end of file diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp index 47ca1fa..640c2d9 100644 --- a/test/testBoundary.cpp +++ b/test/testBoundary.cpp @@ -1,68 +1,74 @@ -#include #include -#include -#include -#include #include +#include +#include +#include +#include -TEST_CASE("BoundaryElement"){ +TEST_CASE("BoundaryElement") { - SUBCASE("Closed case"){ - BoundaryElement boundaryElementClosed = BoundaryElement(); - CHECK_NOTHROW(BoundaryElement()); - CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); - CHECK_EQ(isnan(boundaryElementClosed.getValue()), isnan(NAN)); - CHECK_THROWS(boundaryElementClosed.setValue(0.2)); - } + SUBCASE("Closed case") { + BoundaryElement boundaryElementClosed = BoundaryElement(); + CHECK_NOTHROW(BoundaryElement()); + CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); + CHECK_EQ(isnan(boundaryElementClosed.getValue()), isnan(NAN)); + CHECK_THROWS(boundaryElementClosed.setValue(0.2)); + } - SUBCASE("Constant case"){ - BoundaryElement boundaryElementConstant = BoundaryElement(0.1); - CHECK_NOTHROW(BoundaryElement(0.1)); - CHECK_EQ(boundaryElementConstant.getType(), BC_TYPE_CONSTANT); - CHECK_EQ(boundaryElementConstant.getValue(), 0.1); - CHECK_NOTHROW(boundaryElementConstant.setValue(0.2)); - CHECK_EQ(boundaryElementConstant.getValue(), 0.2); - } + SUBCASE("Constant case") { + BoundaryElement boundaryElementConstant = BoundaryElement(0.1); + CHECK_NOTHROW(BoundaryElement(0.1)); + CHECK_EQ(boundaryElementConstant.getType(), BC_TYPE_CONSTANT); + CHECK_EQ(boundaryElementConstant.getValue(), 0.1); + CHECK_NOTHROW(boundaryElementConstant.setValue(0.2)); + CHECK_EQ(boundaryElementConstant.getValue(), 0.2); + } } -TEST_CASE("Boundary Class"){ - Grid grid1D = Grid(10); - Grid grid2D = Grid(10, 12); - Boundary boundary1D = Boundary(grid1D); - Boundary boundary2D = Boundary(grid2D); - vector boundary1DVector(1, BoundaryElement(1.0)); +TEST_CASE("Boundary Class") { + Grid grid1D = Grid(10); + Grid grid2D = Grid(10, 12); + Boundary boundary1D = Boundary(grid1D); + Boundary boundary2D = Boundary(grid2D); + vector boundary1DVector(1, BoundaryElement(1.0)); - SUBCASE("Boundaries 1D case"){ - CHECK_NOTHROW(Boundary boundary(grid1D)); - CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1); - CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1); - CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CLOSED); - CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_TOP)); - CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_BOTTOM)); - CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT)); - CHECK_THROWS(boundary1D.setBoundarySideClosed(BC_SIDE_TOP)); - CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); - CHECK_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); - CHECK_THROWS(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2)); - CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); - CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); - } + SUBCASE("Boundaries 1D case") { + CHECK_NOTHROW(Boundary boundary(grid1D)); + CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1); + CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1); + CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), + BC_TYPE_CLOSED); + CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_TOP)); + CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_BOTTOM)); + CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT)); + CHECK_THROWS(boundary1D.setBoundarySideClosed(BC_SIDE_TOP)); + CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); + CHECK_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); + CHECK_THROWS(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2)); + CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), + BC_TYPE_CONSTANT); + CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), + boundary1DVector[0].getType()); + } - SUBCASE("Boundaries 2D case"){ - CHECK_NOTHROW(Boundary boundary(grid1D)); - CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10); - CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).size(), 10); - CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_TOP).size(), 12); - CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_BOTTOM).size(), 12); - CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CLOSED); - CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_TOP)); - CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_BOTTOM)); - CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_LEFT)); - CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP)); - CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); - CHECK_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); - CHECK_THROWS(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12)); - CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); - CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); - } + SUBCASE("Boundaries 2D case") { + CHECK_NOTHROW(Boundary boundary(grid1D)); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).size(), 10); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_TOP).size(), 12); + CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_BOTTOM).size(), 12); + CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), + BC_TYPE_CLOSED); + CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_TOP)); + CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_BOTTOM)); + CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_LEFT)); + CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP)); + CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); + CHECK_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); + CHECK_THROWS(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12)); + CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), + BC_TYPE_CONSTANT); + CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), + boundary1DVector[0].getType()); + } } \ No newline at end of file diff --git a/test/testFTCS.cpp b/test/testFTCS.cpp index d77999a..f2cbda2 100644 --- a/test/testFTCS.cpp +++ b/test/testFTCS.cpp @@ -1,19 +1,19 @@ -#include #include <../src/FTCS.cpp> +#include #include TEST_CASE("Maths") { - SUBCASE("mean between two alphas") { - double alpha1 = 10; - double alpha2 = 20; - double average = 15; - double harmonicMean = double(2) / ((double(1)/alpha1)+(double(1)/alpha2)); - - // double difference = std::fabs(calcAlphaIntercell(alpha1, alpha2) - harmonicMean); - // CHECK(difference < std::numeric_limits::epsilon()); - CHECK_EQ(calcAlphaIntercell(alpha1, alpha2), harmonicMean); - CHECK_EQ(calcAlphaIntercell(alpha1, alpha2, false), average); - } - + SUBCASE("mean between two alphas") { + double alpha1 = 10; + double alpha2 = 20; + double average = 15; + double harmonicMean = + double(2) / ((double(1) / alpha1) + (double(1) / alpha2)); + // double difference = std::fabs(calcAlphaIntercell(alpha1, alpha2) - + // harmonicMean); CHECK(difference < + // std::numeric_limits::epsilon()); + CHECK_EQ(calcAlphaIntercell(alpha1, alpha2), harmonicMean); + CHECK_EQ(calcAlphaIntercell(alpha1, alpha2, false), average); + } } \ No newline at end of file diff --git a/test/testGrid.cpp b/test/testGrid.cpp index 4963976..ea16f6f 100644 --- a/test/testGrid.cpp +++ b/test/testGrid.cpp @@ -1,251 +1,249 @@ +#include #include #include -#include TEST_CASE("1D Grid, too small length") { - int l = 2; - CHECK_THROWS(Grid(l)); + int l = 2; + CHECK_THROWS(Grid(l)); } TEST_CASE("2D Grid, too small side") { - int r = 2; - int c = 4; - CHECK_THROWS(Grid(r, c)); + int r = 2; + int c = 4; + CHECK_THROWS(Grid(r, c)); - r = 4; - c = 2; - CHECK_THROWS(Grid(r, c)); + r = 4; + c = 2; + CHECK_THROWS(Grid(r, c)); } TEST_CASE("1D Grid") { - int l = 12; - Grid grid(l); + int l = 12; + Grid grid(l); - SUBCASE("correct construction") { - CHECK_EQ(grid.getDim(), 1); - CHECK_EQ(grid.getLength(), l); - CHECK_EQ(grid.getCol(), l); - CHECK_EQ(grid.getRow(), 1); + SUBCASE("correct construction") { + CHECK_EQ(grid.getDim(), 1); + CHECK_EQ(grid.getLength(), l); + CHECK_EQ(grid.getCol(), l); + CHECK_EQ(grid.getRow(), 1); - CHECK_EQ(grid.getConcentrations().rows(), 1); - CHECK_EQ(grid.getConcentrations().cols(), l); - CHECK_EQ(grid.getAlpha().rows(), 1); - CHECK_EQ(grid.getAlpha().cols(), l); - CHECK_EQ(grid.getDeltaCol(), 1); + CHECK_EQ(grid.getConcentrations().rows(), 1); + CHECK_EQ(grid.getConcentrations().cols(), l); + CHECK_EQ(grid.getAlpha().rows(), 1); + CHECK_EQ(grid.getAlpha().cols(), l); + CHECK_EQ(grid.getDeltaCol(), 1); - CHECK_THROWS(grid.getAlphaX()); - CHECK_THROWS(grid.getAlphaY()); - CHECK_THROWS(grid.getDeltaRow()); - } + CHECK_THROWS(grid.getAlphaX()); + CHECK_THROWS(grid.getAlphaY()); + CHECK_THROWS(grid.getDeltaRow()); + } - SUBCASE("setting concentrations") { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(1, l, 3); - CHECK_NOTHROW(grid.setConcentrations(concentrations)); + SUBCASE("setting concentrations") { + // correct concentrations matrix + MatrixXd concentrations = MatrixXd::Constant(1, l, 3); + CHECK_NOTHROW(grid.setConcentrations(concentrations)); - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - } + // false concentrations matrix + MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + } - SUBCASE("setting alpha") { - // correct alpha matrix - MatrixXd alpha = MatrixXd::Constant(1, l, 3); - CHECK_NOTHROW(grid.setAlpha(alpha)); + SUBCASE("setting alpha") { + // correct alpha matrix + MatrixXd alpha = MatrixXd::Constant(1, l, 3); + CHECK_NOTHROW(grid.setAlpha(alpha)); - CHECK_THROWS(grid.setAlpha(alpha, alpha)); + CHECK_THROWS(grid.setAlpha(alpha, alpha)); - grid.setAlpha(alpha); - CHECK_EQ(grid.getAlpha(), alpha); - CHECK_THROWS(grid.getAlphaX()); - CHECK_THROWS(grid.getAlphaY()); + grid.setAlpha(alpha); + CHECK_EQ(grid.getAlpha(), alpha); + CHECK_THROWS(grid.getAlphaX()); + CHECK_THROWS(grid.getAlphaY()); - // false alpha matrix - MatrixXd wAlpha = MatrixXd::Constant(3, l, 2); - CHECK_THROWS(grid.setAlpha(wAlpha)); - } + // false alpha matrix + MatrixXd wAlpha = MatrixXd::Constant(3, l, 2); + CHECK_THROWS(grid.setAlpha(wAlpha)); + } - SUBCASE("setting domain") { - int d = 8; - // set 1D domain - CHECK_NOTHROW(grid.setDomain(d)); + SUBCASE("setting domain") { + int d = 8; + // set 1D domain + CHECK_NOTHROW(grid.setDomain(d)); - // set 2D domain - CHECK_THROWS(grid.setDomain(d, d)); + // set 2D domain + CHECK_THROWS(grid.setDomain(d, d)); - grid.setDomain(d); - CHECK_EQ(grid.getDeltaCol(), double(d)/double(l)); - CHECK_THROWS(grid.getDeltaRow()); + grid.setDomain(d); + CHECK_EQ(grid.getDeltaCol(), double(d) / double(l)); + CHECK_THROWS(grid.getDeltaRow()); - // set too small domain - d = 0; - CHECK_THROWS(grid.setDomain(d)); - d = -2; - CHECK_THROWS(grid.setDomain(d)); - } + // set too small domain + d = 0; + CHECK_THROWS(grid.setDomain(d)); + d = -2; + CHECK_THROWS(grid.setDomain(d)); + } } TEST_CASE("2D Grid quadratic") { - int rc = 12; - Grid grid(rc, rc); + int rc = 12; + Grid grid(rc, rc); - SUBCASE("correct construction") { - CHECK_EQ(grid.getDim(), 2); - CHECK_THROWS(grid.getLength()); - CHECK_EQ(grid.getCol(), rc); - CHECK_EQ(grid.getRow(), rc); + SUBCASE("correct construction") { + CHECK_EQ(grid.getDim(), 2); + CHECK_THROWS(grid.getLength()); + CHECK_EQ(grid.getCol(), rc); + CHECK_EQ(grid.getRow(), rc); - CHECK_EQ(grid.getConcentrations().rows(), rc); - CHECK_EQ(grid.getConcentrations().cols(), rc); - CHECK_THROWS(grid.getAlpha()); + CHECK_EQ(grid.getConcentrations().rows(), rc); + CHECK_EQ(grid.getConcentrations().cols(), rc); + CHECK_THROWS(grid.getAlpha()); - CHECK_EQ(grid.getAlphaX().rows(), rc); - CHECK_EQ(grid.getAlphaX().cols(), rc); - CHECK_EQ(grid.getAlphaY().rows(), rc); - CHECK_EQ(grid.getAlphaY().cols(), rc); - CHECK_EQ(grid.getDeltaRow(), 1); - CHECK_EQ(grid.getDeltaCol(), 1); - } + CHECK_EQ(grid.getAlphaX().rows(), rc); + CHECK_EQ(grid.getAlphaX().cols(), rc); + CHECK_EQ(grid.getAlphaY().rows(), rc); + CHECK_EQ(grid.getAlphaY().cols(), rc); + CHECK_EQ(grid.getDeltaRow(), 1); + CHECK_EQ(grid.getDeltaCol(), 1); + } - SUBCASE("setting concentrations") { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2); - CHECK_NOTHROW(grid.setConcentrations(concentrations)); + SUBCASE("setting concentrations") { + // correct concentrations matrix + MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2); + CHECK_NOTHROW(grid.setConcentrations(concentrations)); + // false concentrations matrix + MatrixXd wConcentrations = MatrixXd::Constant(rc, rc + 3, 1); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + wConcentrations = MatrixXd::Constant(rc + 3, rc, 4); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + wConcentrations = MatrixXd::Constant(rc + 2, rc + 4, 6); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + } - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(rc, rc+3, 1); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(rc+3, rc, 4); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(rc+2, rc+4, 6); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - } + SUBCASE("setting alphas") { + // correct alpha matrices + MatrixXd alphax = MatrixXd::Constant(rc, rc, 2); + MatrixXd alphay = MatrixXd::Constant(rc, rc, 4); + CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); - SUBCASE("setting alphas") { - // correct alpha matrices - MatrixXd alphax = MatrixXd::Constant(rc, rc, 2); - MatrixXd alphay = MatrixXd::Constant(rc, rc, 4); - CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); + CHECK_THROWS(grid.setAlpha(alphax)); - CHECK_THROWS(grid.setAlpha(alphax)); + grid.setAlpha(alphax, alphay); + CHECK_EQ(grid.getAlphaX(), alphax); + CHECK_EQ(grid.getAlphaY(), alphay); + CHECK_THROWS(grid.getAlpha()); - grid.setAlpha(alphax, alphay); - CHECK_EQ(grid.getAlphaX(), alphax); - CHECK_EQ(grid.getAlphaY(), alphay); - CHECK_THROWS(grid.getAlpha()); + // false alpha matrices + alphax = MatrixXd::Constant(rc + 3, rc + 1, 3); + CHECK_THROWS(grid.setAlpha(alphax, alphay)); + alphay = MatrixXd::Constant(rc + 2, rc + 1, 3); + CHECK_THROWS(grid.setAlpha(alphax, alphay)); + } - // false alpha matrices - alphax = MatrixXd::Constant(rc+3, rc+1, 3); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); - alphay = MatrixXd::Constant(rc+2, rc+1, 3); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); - } + SUBCASE("setting domain") { + int dr = 8; + int dc = 9; - SUBCASE("setting domain") { - int dr = 8; - int dc = 9; + // set 1D domain + CHECK_THROWS(grid.setDomain(dr)); - // set 1D domain - CHECK_THROWS(grid.setDomain(dr)); + // set 2D domain + CHECK_NOTHROW(grid.setDomain(dr, dc)); - // set 2D domain - CHECK_NOTHROW(grid.setDomain(dr, dc)); + grid.setDomain(dr, dc); + CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc)); + CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc)); - grid.setDomain(dr, dc); - CHECK_EQ(grid.getDeltaCol(), double(dc)/double(rc)); - CHECK_EQ(grid.getDeltaRow(), double(dr)/double(rc)); - - // set too small domain - dr = 0; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = 8; - dc = 0; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = -2; - CHECK_THROWS(grid.setDomain(dr, dc)); - } + // set too small domain + dr = 0; + CHECK_THROWS(grid.setDomain(dr, dc)); + dr = 8; + dc = 0; + CHECK_THROWS(grid.setDomain(dr, dc)); + dr = -2; + CHECK_THROWS(grid.setDomain(dr, dc)); + } } TEST_CASE("2D Grid non-quadratic") { - int r = 12; - int c = 15; - Grid grid(r, c); + int r = 12; + int c = 15; + Grid grid(r, c); - SUBCASE("correct construction") { - CHECK_EQ(grid.getDim(), 2); - CHECK_THROWS(grid.getLength()); - CHECK_EQ(grid.getCol(), c); - CHECK_EQ(grid.getRow(), r); + SUBCASE("correct construction") { + CHECK_EQ(grid.getDim(), 2); + CHECK_THROWS(grid.getLength()); + CHECK_EQ(grid.getCol(), c); + CHECK_EQ(grid.getRow(), r); - CHECK_EQ(grid.getConcentrations().rows(), r); - CHECK_EQ(grid.getConcentrations().cols(), c); - CHECK_THROWS(grid.getAlpha()); + CHECK_EQ(grid.getConcentrations().rows(), r); + CHECK_EQ(grid.getConcentrations().cols(), c); + CHECK_THROWS(grid.getAlpha()); - CHECK_EQ(grid.getAlphaX().rows(), r); - CHECK_EQ(grid.getAlphaX().cols(), c); - CHECK_EQ(grid.getAlphaY().rows(), r); - CHECK_EQ(grid.getAlphaY().cols(), c); - CHECK_EQ(grid.getDeltaRow(), 1); - CHECK_EQ(grid.getDeltaCol(), 1); - } + CHECK_EQ(grid.getAlphaX().rows(), r); + CHECK_EQ(grid.getAlphaX().cols(), c); + CHECK_EQ(grid.getAlphaY().rows(), r); + CHECK_EQ(grid.getAlphaY().cols(), c); + CHECK_EQ(grid.getDeltaRow(), 1); + CHECK_EQ(grid.getDeltaCol(), 1); + } - SUBCASE("setting concentrations") { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(r, c, 2); - CHECK_NOTHROW(grid.setConcentrations(concentrations)); + SUBCASE("setting concentrations") { + // correct concentrations matrix + MatrixXd concentrations = MatrixXd::Constant(r, c, 2); + CHECK_NOTHROW(grid.setConcentrations(concentrations)); + // false concentrations matrix + MatrixXd wConcentrations = MatrixXd::Constant(r, c + 3, 6); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + wConcentrations = MatrixXd::Constant(r + 3, c, 3); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + wConcentrations = MatrixXd::Constant(r + 2, c + 4, 2); + CHECK_THROWS(grid.setConcentrations(wConcentrations)); + } - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(r, c+3, 6); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(r+3, c, 3); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(r+2, c+4, 2); - CHECK_THROWS(grid.setConcentrations(wConcentrations)); - } + SUBCASE("setting alphas") { + // correct alpha matrices + MatrixXd alphax = MatrixXd::Constant(r, c, 2); + MatrixXd alphay = MatrixXd::Constant(r, c, 4); + CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); - SUBCASE("setting alphas") { - // correct alpha matrices - MatrixXd alphax = MatrixXd::Constant(r, c, 2); - MatrixXd alphay = MatrixXd::Constant(r, c, 4); - CHECK_NOTHROW(grid.setAlpha(alphax, alphay)); + CHECK_THROWS(grid.setAlpha(alphax)); - CHECK_THROWS(grid.setAlpha(alphax)); + grid.setAlpha(alphax, alphay); + CHECK_EQ(grid.getAlphaX(), alphax); + CHECK_EQ(grid.getAlphaY(), alphay); + CHECK_THROWS(grid.getAlpha()); - grid.setAlpha(alphax, alphay); - CHECK_EQ(grid.getAlphaX(), alphax); - CHECK_EQ(grid.getAlphaY(), alphay); - CHECK_THROWS(grid.getAlpha()); + // false alpha matrices + alphax = MatrixXd::Constant(r + 3, c + 1, 3); + CHECK_THROWS(grid.setAlpha(alphax, alphay)); + alphay = MatrixXd::Constant(r + 2, c + 1, 5); + CHECK_THROWS(grid.setAlpha(alphax, alphay)); + } - // false alpha matrices - alphax = MatrixXd::Constant(r+3, c+1, 3); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); - alphay = MatrixXd::Constant(r+2, c+1, 5); - CHECK_THROWS(grid.setAlpha(alphax, alphay)); - } + SUBCASE("setting domain") { + int dr = 8; + int dc = 9; - SUBCASE("setting domain") { - int dr = 8; - int dc = 9; + // set 1D domain + CHECK_THROWS(grid.setDomain(dr)); - // set 1D domain - CHECK_THROWS(grid.setDomain(dr)); + // set 2D domain + CHECK_NOTHROW(grid.setDomain(dr, dc)); - // set 2D domain - CHECK_NOTHROW(grid.setDomain(dr, dc)); + grid.setDomain(dr, dc); + CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c)); + CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r)); - grid.setDomain(dr, dc); - CHECK_EQ(grid.getDeltaCol(), double(dc)/double(c)); - CHECK_EQ(grid.getDeltaRow(), double(dr)/double(r)); - - // set too small domain - dr = 0; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = 8; - dc = -1; - CHECK_THROWS(grid.setDomain(dr, dc)); - dr = -2; - CHECK_THROWS(grid.setDomain(dr, dc)); - } + // set too small domain + dr = 0; + CHECK_THROWS(grid.setDomain(dr, dc)); + dr = 8; + dc = -1; + CHECK_THROWS(grid.setDomain(dr, dc)); + dr = -2; + CHECK_THROWS(grid.setDomain(dr, dc)); + } } \ No newline at end of file diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp index 4efc1b0..6a661e4 100644 --- a/test/testSimulation.cpp +++ b/test/testSimulation.cpp @@ -1,110 +1,104 @@ -#include -#include -#include #include "TestUtils.cpp" +#include +#include #include +#include // include the configured header file #include -static Grid setupSimulation(APPROACH approach, double timestep, int iterations) { - int row = 11; - int col = 11; - int domain_row = 10; - int domain_col = 10; +static Grid setupSimulation(APPROACH approach, double timestep, + int iterations) { + int row = 11; + int col = 11; + int domain_row = 10; + int domain_col = 10; + // Grid + Grid grid = Grid(row, col); + grid.setDomain(domain_row, domain_col); - // Grid - Grid grid = Grid(row, col); - grid.setDomain(domain_row, domain_col); + MatrixXd concentrations = MatrixXd::Constant(row, col, 0); + concentrations(5, 5) = 1; + grid.setConcentrations(concentrations); - MatrixXd concentrations = MatrixXd::Constant(row, col, 0); - concentrations(5,5) = 1; - grid.setConcentrations(concentrations); - - MatrixXd alpha = MatrixXd::Constant(row, col, 1); - for (int i = 0; i < 5; i++) { - for (int j = 0; j < 6; j++) { - alpha(i, j) = 0.01; - } + MatrixXd alpha = MatrixXd::Constant(row, col, 1); + for (int i = 0; i < 5; i++) { + for (int j = 0; j < 6; j++) { + alpha(i, j) = 0.01; } - for (int i = 0; i < 5; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.001; - } + } + for (int i = 0; i < 5; i++) { + for (int j = 6; j < 11; j++) { + alpha(i, j) = 0.001; } - for (int i = 5; i < 11; i++) { - for (int j = 6; j < 11; j++) { - alpha(i, j) = 0.1; - } + } + for (int i = 5; i < 11; i++) { + for (int j = 6; j < 11; j++) { + alpha(i, j) = 0.1; } - grid.setAlpha(alpha, alpha); + } + grid.setAlpha(alpha, alpha); + // Boundary + Boundary bc = Boundary(grid); - // Boundary - Boundary bc = Boundary(grid); + // Simulation + Simulation sim = Simulation(grid, bc, approach); + // sim.setOutputConsole(CONSOLE_OUTPUT_ON); + sim.setTimestep(timestep); + sim.setIterations(iterations); + sim.run(); - - // Simulation - Simulation sim = Simulation(grid, bc, approach); - // sim.setOutputConsole(CONSOLE_OUTPUT_ON); - sim.setTimestep(timestep); - sim.setIterations(iterations); - sim.run(); - - // RUN - return grid; - + // RUN + return grid; } TEST_CASE("equality to reference matrix with FTCS") { - // set string from the header file - string test_path = testSimulationCSVDir; - MatrixXd reference = CSV2Eigen(test_path); - cout << "FTCS Test: " << endl; - Grid grid = setupSimulation(FTCS_APPROACH, 0.001, 7000); - cout << endl; - CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true); + // set string from the header file + string test_path = testSimulationCSVDir; + MatrixXd reference = CSV2Eigen(test_path); + cout << "FTCS Test: " << endl; + Grid grid = setupSimulation(FTCS_APPROACH, 0.001, 7000); + cout << endl; + 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); - cout << "BTCS Test: " << endl; - Grid grid = setupSimulation(BTCS_APPROACH, 1, 7); - cout << endl; - CHECK(checkSimilarityV2(reference, grid.getConcentrations(), 0.01) == true); + // set string from the header file + string test_path = testSimulationCSVDir; + MatrixXd reference = CSV2Eigen(test_path); + cout << "BTCS Test: " << endl; + Grid grid = setupSimulation(BTCS_APPROACH, 1, 7); + cout << endl; + CHECK(checkSimilarityV2(reference, grid.getConcentrations(), 0.01) == true); } -TEST_CASE("Initialize environment"){ - int rc = 12; - Grid grid(rc, rc); - Boundary boundary(grid); +TEST_CASE("Initialize environment") { + int rc = 12; + Grid grid(rc, rc); + Boundary boundary(grid); - CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH)); + CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH)); } -TEST_CASE("Simulation environment"){ - int rc = 12; - Grid grid(rc, rc); - Boundary boundary(grid); - Simulation sim(grid, boundary, FTCS_APPROACH); +TEST_CASE("Simulation environment") { + int rc = 12; + Grid grid(rc, rc); + Boundary boundary(grid); + Simulation sim(grid, boundary, FTCS_APPROACH); - SUBCASE("default paremeters") { - CHECK_EQ(sim.getIterations(), -1); - } + SUBCASE("default paremeters") { CHECK_EQ(sim.getIterations(), -1); } - SUBCASE("set iterations") { - CHECK_NOTHROW(sim.setIterations(2000)); - CHECK_EQ(sim.getIterations(), 2000); - CHECK_THROWS(sim.setIterations(-300)); - } + SUBCASE("set iterations") { + CHECK_NOTHROW(sim.setIterations(2000)); + CHECK_EQ(sim.getIterations(), 2000); + CHECK_THROWS(sim.setIterations(-300)); + } - SUBCASE("set timestep") { - CHECK_NOTHROW(sim.setTimestep(0.1)); - CHECK_EQ(sim.getTimestep(), 0.1); - CHECK_THROWS(sim.setTimestep(-0.3)); - } + SUBCASE("set timestep") { + CHECK_NOTHROW(sim.setTimestep(0.1)); + CHECK_EQ(sim.getTimestep(), 0.1); + CHECK_THROWS(sim.setTimestep(-0.3)); + } } -