From 77914ea69faffe8cddcdcfd6a8e7c7172566d748 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 17 Oct 2023 12:07:33 +0200 Subject: [PATCH] fix: include optional output of csv during thomas algorithm fix: marco's benchmark --- include/tug/Core/BTCS.hpp | 23 +++++ naaice/CMakeLists.txt | 2 + naaice/NAAICE_dble_vs_float.cpp | 168 +++++++++++--------------------- 3 files changed, 82 insertions(+), 111 deletions(-) diff --git a/include/tug/Core/BTCS.hpp b/include/tug/Core/BTCS.hpp index cea830a..b844ba4 100644 --- a/include/tug/Core/BTCS.hpp +++ b/include/tug/Core/BTCS.hpp @@ -280,6 +280,29 @@ static Eigen::VectorX ThomasAlgorithm(Eigen::SparseMatrix &A, a_diag[n - 1] = A.coeff(n - 1, n - 2); b_diag[n - 1] = A.coeff(n - 1, n - 1); + // HACK: write CSV to file +#ifdef WRITE_THOMAS_CSV +#include +#include + static std::uint32_t file_index = 0; + std::string file_name = "Thomas_" + std::to_string(file_index++) + ".csv"; + + std::ofstream out_file; + + out_file.open(file_name, std::ofstream::trunc | std::ofstream::out); + + // print header + out_file << "Aa, Ab, Ac, b\n"; + + // iterate through all elements + for (std::size_t i = 0; i < n; i++) { + out_file << a_diag[i] << ", " << b_diag[i] << ", " << c_diag[i] << ", " + << b[i] << "\n"; + } + + out_file.close(); +#endif + // start solving - c_diag and x_vec are overwritten n--; c_diag[0] /= b_diag[0]; diff --git a/naaice/CMakeLists.txt b/naaice/CMakeLists.txt index b0a45d0..209dc91 100644 --- a/naaice/CMakeLists.txt +++ b/naaice/CMakeLists.txt @@ -8,5 +8,7 @@ configure_file(files.hpp.in files.hpp) target_include_directories(naaice PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) target_include_directories(NAAICE_dble_vs_float PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) +target_compile_definitions(naaice PRIVATE WRITE_THOMAS_CSV) + target_link_libraries(naaice PUBLIC tug) target_link_libraries(NAAICE_dble_vs_float PUBLIC tug) diff --git a/naaice/NAAICE_dble_vs_float.cpp b/naaice/NAAICE_dble_vs_float.cpp index 5be74a4..ac3718c 100644 --- a/naaice/NAAICE_dble_vs_float.cpp +++ b/naaice/NAAICE_dble_vs_float.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -7,8 +8,8 @@ #include #include #include +#include #include -#include #include "files.hpp" @@ -96,151 +97,96 @@ Eigen::MatrixXd rmVecTocmMatrix(const std::vector> &vec, return out_mat; } +template int doWork(int ngrid) { + + constexpr T dt = 10; -int DoDble(int ngrid, APPROACH approach) { - - constexpr double dt = 10; - // create a grid - std::cout << "DOUBLE grid: " << ngrid << std::endl; + std::string name; - Grid64 grid64(ngrid, ngrid); + if constexpr (std::is_same_v) { + name = "DOUBLE"; + } else if constexpr (std::is_same_v) { + name = "FLOAT"; + } else { + name = "unknown"; + } + + std::cout << name << " grid: " << ngrid << std::endl; + + Grid grid(ngrid, ngrid); + // Grid64 grid(ngrid, ngrid); // (optional) set the domain, e.g.: - grid64.setDomain(0.1, 0.1); + grid.setDomain(0.1, 0.1); - Eigen::MatrixXd initConc64 = Eigen::MatrixXd::Constant(ngrid, ngrid, 0); + Eigen::MatrixX initConc64 = Eigen::MatrixX::Constant(ngrid, ngrid, 0); initConc64(50, 50) = 1E-6; - grid64.setConcentrations(initConc64); - - const double sum_init64 = initConc64.sum(); + grid.setConcentrations(initConc64); - constexpr double alphax_val = 5e-10; - Eigen::MatrixXd alphax = Eigen::MatrixXd::Constant(ngrid, ngrid, alphax_val); // row,col,value - constexpr double alphay_val = 1e-10; - Eigen::MatrixXd alphay = Eigen::MatrixXd::Constant(ngrid, ngrid, alphay_val); // row,col,value - grid64.setAlpha(alphax, alphay); + const T sum_init64 = initConc64.sum(); + + constexpr T alphax_val = 5e-10; + Eigen::MatrixX alphax = + Eigen::MatrixX::Constant(ngrid, ngrid, alphax_val); // row,col,value + constexpr T alphay_val = 1e-10; + Eigen::MatrixX alphay = + Eigen::MatrixX::Constant(ngrid, ngrid, alphay_val); // row,col,value + grid.setAlpha(alphax, alphay); // create a boundary with constant values - Boundary bc64 = Boundary(grid64); - bc64.setBoundarySideClosed(BC_SIDE_LEFT); - bc64.setBoundarySideClosed(BC_SIDE_RIGHT); - bc64.setBoundarySideClosed(BC_SIDE_TOP); - bc64.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); // set up a simulation environment - Simulation Sim64 = - Simulation(grid64, bc64, approach); // grid_64,boundary,simulation-approach + Simulation Sim = + Simulation(grid, bc); // grid_64,boundary,simulation-approach + + // Sim64.setSolver(THOMAS_ALGORITHM_SOLVER); - //Sim64.setSolver(THOMAS_ALGORITHM_SOLVER); - // set the timestep of the simulation - Sim64.setTimestep(dt); // timestep + Sim.setTimestep(dt); // timestep // set the number of iterations - Sim64.setIterations(2); + Sim.setIterations(2); // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, // CSV_OUTPUT_VERBOSE] - Sim64.setOutputCSV(CSV_OUTPUT_ON); + Sim.setOutputCSV(CSV_OUTPUT_ON); // set output to the console to 'ON' - Sim64.setOutputConsole(CONSOLE_OUTPUT_OFF); + Sim.setOutputConsole(CONSOLE_OUTPUT_OFF); // // **** RUN SIM64 **** - auto begin64 = std::chrono::high_resolution_clock::now(); + auto begin_t = std::chrono::high_resolution_clock::now(); - Sim64.run(); + Sim.run(); - auto end64 = std::chrono::high_resolution_clock::now(); - auto ms64 = std::chrono::duration_cast(end64 - begin64); + auto end_t = std::chrono::high_resolution_clock::now(); + auto ms_t = + std::chrono::duration_cast(end_t - begin_t); - const double sum_after64 = grid64.getConcentrations().sum(); + const double sum_after64 = grid.getConcentrations().sum(); - std::cout << "Sum of init field: " << std::setprecision(15) << sum_init64 - << "\nSum after 2 iterations: " << sum_after64 - << "\nMilliseconds: " << ms64.count() << std::endl << std::endl; + std::cout << "Sum of init field: " << std::setprecision(15) << sum_init64 + << "\nSum after 2 iterations: " << sum_after64 + << "\nMilliseconds: " << ms_t.count() << std::endl + << std::endl; return 0; - } -int DoFloat(int ngrid, APPROACH approach) { - - constexpr double dt = 10; - - // create a grid - std::cout << "FLOAT grid: " << ngrid << std::endl; - - Grid32 grid32(ngrid, ngrid); - - // (optional) set the domain, e.g.: - grid32.setDomain(0.1, 0.1); - - Eigen::MatrixXf initConc32 = Eigen::MatrixXf::Constant(ngrid, ngrid, 0); - initConc32(50, 50) = 1E-6; - grid32.setConcentrations(initConc32); - - const float sum_init32 = initConc32.sum(); - - constexpr float alphax_val = 5e-10; - Eigen::MatrixXf alphax = Eigen::MatrixXf::Constant(ngrid, ngrid, alphax_val); // row,col,value - constexpr float alphay_val = 1e-10; - Eigen::MatrixXf alphay = Eigen::MatrixXf::Constant(ngrid, ngrid, alphay_val); // row,col,value - grid32.setAlpha(alphax, alphay); - - // create a boundary with constant values - Boundary bc32 = Boundary(grid32); - bc32.setBoundarySideClosed(BC_SIDE_LEFT); - bc32.setBoundarySideClosed(BC_SIDE_RIGHT); - bc32.setBoundarySideClosed(BC_SIDE_TOP); - bc32.setBoundarySideClosed(BC_SIDE_BOTTOM); - - // set up a simulation environment - Simulation Sim32 = - Simulation(grid32, bc32, approach); // grid_32,boundary,simulation-approach - - // Sim32.setSolver(THOMAS_ALGORITHM_SOLVER); - - // set the timestep of the simulation - Sim32.setTimestep(dt); // timestep - - // set the number of iterations - Sim32.setIterations(2); - - // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, - // CSV_OUTPUT_VERBOSE] - Sim32.setOutputCSV(CSV_OUTPUT_ON); - - // set output to the console to 'ON' - Sim32.setOutputConsole(CONSOLE_OUTPUT_OFF); - - // // **** RUN SIM32 **** - auto begin32 = std::chrono::high_resolution_clock::now(); - - Sim32.run(); - - auto end32 = std::chrono::high_resolution_clock::now(); - auto ms32 = std::chrono::duration_cast(end32 - begin32); - - const float sum_after32 = grid32.getConcentrations().sum(); - - std::cout << "Sum of init field: " << std::setprecision(15) << sum_init32 - << "\nSum after 2 iteration: " << sum_after32 - << "\nMilliseconds: " << ms32.count() << std::endl << std::endl; - return 0; - -} - - int main(int argc, char *argv[]) { int n[] = {101, 201, 501, 1001, 2001}; for (int i = 0; i < std::size(n); i++) { - DoFloat(n[i], BTCS_APPROACH); - DoDble(n[i], BTCS_APPROACH); - DoFloat(n[i], FTCS_APPROACH); - DoDble(n[i], FTCS_APPROACH); + doWork(n[i]); + doWork(n[i]); + doWork(n[i]); + doWork(n[i]); } return 0;