diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h index 04670af..c773f15 100644 --- a/nufi/nufi_solver.h +++ b/nufi/nufi_solver.h @@ -11,7 +11,6 @@ #include "nufi/fields.h" // dont remove #include "nufi/grids.h" #include "nufi/parameters.h" -#include "nufi/poisson_problem.h" using namespace dealii; @@ -19,7 +18,6 @@ class NuFISolver { public: NuFISolver(); - void run(); std::vector eval_rho(unsigned int n, std::vector &x, const std::vector> &grid_struct, @@ -48,10 +46,6 @@ private: double x_min = Parameters::X_DOMAIN_LEFT; double x_max = Parameters::X_DOMAIN_RIGHT; - std::vector rho; - unsigned int order; - - PoissonProblem<1> poisson; }; #endif diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h index 2a6f261..1cb3713 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -81,7 +81,6 @@ public: unsigned int get_dof_size(); void set_rhs_function( std::function(const std::vector> &)> f); - // void set_rhs(const Vector &new_rhs) { rhs = new_rhs; } const Vector &get_solution() const { return solution; } const MappingQ &get_mapping() const { return mapping; } diff --git a/src/main.cc b/src/main.cc index b6f9e42..8b60171 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,6 +1,8 @@ #include #include #include +#include +#include #include void clear_results_directory(const std::string &dir) { @@ -13,16 +15,132 @@ void clear_results_directory(const std::string &dir) { std::cout << "Deleted: " << entry.path() << '\n'; } } -} +}; + +template void run() { + unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT); + + PoissonProblem poisson(dim); + NuFISolver solver; + using std::abs; + using std::max; + + std::cout << "Initializing dealii Poisson Solver\n"; + poisson.initialize(); + + std::vector int_E_squared; + int_E_squared.reserve(Nt); + std::vector int_E_squared_times; + int_E_squared_times.reserve(Nt); + + std::vector> grid_versions; + std::vector> phi_history; + + std::ofstream time_file(Parameters::PLOT_DIR + "simulation_time.dat"); + + double total_time = 0; + stopwatch total_timer; + + time_file << "it " + << "step_time " + << "total_time " + << "compute_time " + << "refine_time " + << "plot_time" + << "\n"; + + std::ofstream error_file(Parameters::PLOT_DIR + "error_estimate.dat"); + error_file << "# nufi Kelly l2 error estimate\n"; + error_file << "# it l2_error_estimate\n"; + + for (unsigned int it = 0; it < Nt; ++it) { + stopwatch timer; + + double time_elapsed_before = timer.elapsed(); + double compute_time = 0.0; + double refine_time = 0.0; + double plot_time = 0.0; + + std::cout << "Timestep " << it << " / " << Nt + << " (simulation time = " << it * Parameters::DT << ")" + << "\n"; + + // START: diagnostics + std::cout << "cells = " << poisson.get_triangulation().n_active_cells() + << "\n" + << " dofs = " << poisson.get_dof_handler().n_dofs() << "\n"; + // END: diagnostics + + double compute_start = timer.elapsed(); + // compute rho + poisson.set_rhs_function([&](const std::vector> &points) { + std::vector x(points.size()); + + for (size_t i = 0; i < points.size(); ++i) + x[i] = points[i][0]; + + return solver.eval_rho(it, x, grid_versions, phi_history, Parameters::NV); + }); + + if (it % Parameters::REFINE_FREQUENCY == 0) { + refine_time = poisson.solve_step(it, grid_versions, true); + compute_time = timer.elapsed() - compute_start - refine_time; + } else { + poisson.solve_step(it, grid_versions, false); + compute_time = timer.elapsed() - compute_start; + } + update_solution_history(phi_history, poisson, + grid_versions.back().grid_version); + + error_file << it << " " << poisson.get_error_estimate() << "\n"; + error_file.flush(); + + double timer_elapsed = timer.elapsed(); + double step_time = timer_elapsed - time_elapsed_before; + + std::cout << "step made in " << step_time << " seconds\n\n"; + + if (it % Parameters::PLOT_FREQUENCY == 0) { + double plot_start = timer.elapsed(); + + std::cout << "Saving results... "; + DiagnosticsSnapshot snap = + compute_diagnostics(solver, it, grid_versions, phi_history, + Parameters::PLOT_NX, Parameters::NV); + + save_f(snap, Parameters::PLOT_DIR + "f_" + std::to_string(it) + ".dat"); + save_rho(snap, + Parameters::PLOT_DIR + "rho_" + std::to_string(it) + ".dat"); + save_Efield(snap, + Parameters::PLOT_DIR + "E_" + std::to_string(it) + ".dat"); + + int_E_squared.push_back(compute_int_E_squared(snap)); + int_E_squared_times.push_back(it * Parameters::DT); + save_time_series(int_E_squared_times, int_E_squared, + Parameters::PLOT_DIR + "int_E_sqr.dat"); + + plot_time = timer.elapsed() - plot_start; + std::cout << "Results saved in " << plot_start << "[s]" << "\n"; + } + + total_time = total_timer.elapsed(); + std::cout << "Time since start = " << total_time << "\n\n"; + + time_file << it << " " << step_time << " " << total_time << " " + << compute_time << " " << refine_time << " " << plot_time << "\n"; + time_file.flush(); + } + + std::cout << "NuFI simulation finished in " << total_time << " seconds.\n"; +}; int main() { // omp_set_max_active_levels(1); std::cout << "Threads: " << omp_get_max_threads() << "\n"; try { - clear_results_directory("results"); - NuFISolver solver; - solver.run(); + clear_results_directory("results"); + run<1>(); // 1 = space_dim } catch (const std::exception &exc) { std::cerr << "\nException:\n" << exc.what() << "\n"; diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index 2a352d5..480f328 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -183,129 +183,4 @@ NuFISolver::eval_rho_points(unsigned int n, const std::vector> &points, return NuFISolver::eval_rho(n, point_vector, grid_struct, phi_history, Nv); } -void NuFISolver::run() { - - //====//====// - // Run prep // - //====//====// - - using std::abs; - using std::max; - - std::vector int_E_squared; - int_E_squared.reserve(Nt); - std::vector int_E_squared_times; - int_E_squared_times.reserve(Nt); - - std::vector> grid_versions; - std::vector> phi_history; - - std::ofstream time_file(Parameters::PLOT_DIR + "simulation_time.dat"); - - double total_time = 0; - stopwatch total_timer; - - time_file << "it " - << "step_time " - << "total_time " - << "compute_time " - << "refine_time " - << "plot_time" - << "\n"; - - std::ofstream error_file(Parameters::PLOT_DIR + "error_estimate.dat"); - error_file << "# nufi Kelly l2 error estimate\n"; - error_file << "# it l2_error_estimate\n"; - - //====//====// - // Time loop// - //====//====// - for (unsigned int it = 0; it < Nt; ++it) { - stopwatch timer; - - double time_elapsed_before = timer.elapsed(); - double compute_time = 0.0; - double refine_time = 0.0; - double plot_time = 0.0; - - std::cout << "Timestep " << it << " / " << Nt - << " (simulation time = " << it * Parameters::DT << ")" - << "\n"; - - // START: diagnostics - std::cout << "cells = " << poisson.get_triangulation().n_active_cells() - << "\n" - << " dofs = " << poisson.get_dof_handler().n_dofs() << "\n"; - // END: diagnostics - - double compute_start = timer.elapsed(); - // compute rho - poisson.set_rhs_function([&](const std::vector> &points) { - std::vector x(points.size()); - - for (size_t i = 0; i < points.size(); ++i) - x[i] = points[i][0]; - - return eval_rho(it, x, grid_versions, phi_history, Parameters::NV); - }); - - if (it % Parameters::REFINE_FREQUENCY == 0) { - refine_time = poisson.solve_step(it, grid_versions, true); - compute_time = timer.elapsed() - compute_start - refine_time; - } else { - poisson.solve_step(it, grid_versions, false); - compute_time = timer.elapsed() - compute_start; - } - update_solution_history(phi_history, poisson, - grid_versions.back().grid_version); - - error_file << it << " " << poisson.get_error_estimate() << "\n"; - error_file.flush(); - - double timer_elapsed = timer.elapsed(); - double step_time = timer_elapsed - time_elapsed_before; - - std::cout << "step made in " << step_time << " seconds\n\n"; - - //====//====// - // Plotting // - //====//====// - - if (it % Parameters::PLOT_FREQUENCY == 0) { - double plot_start = timer.elapsed(); - - std::cout << "Saving results... "; - DiagnosticsSnapshot snap = - compute_diagnostics(*this, it, grid_versions, phi_history, - Parameters::PLOT_NX, Parameters::NV); - - save_f(snap, Parameters::PLOT_DIR + "f_" + std::to_string(it) + ".dat"); - save_rho(snap, - Parameters::PLOT_DIR + "rho_" + std::to_string(it) + ".dat"); - save_Efield(snap, - Parameters::PLOT_DIR + "E_" + std::to_string(it) + ".dat"); - - int_E_squared.push_back(compute_int_E_squared(snap)); - int_E_squared_times.push_back(it * Parameters::DT); - save_time_series(int_E_squared_times, int_E_squared, - Parameters::PLOT_DIR + "int_E_sqr.dat"); - - plot_time = timer.elapsed() - plot_start; - std::cout << "Results saved in " << plot_start << "[s]" << "\n"; - } - - total_time = total_timer.elapsed(); - std::cout << "Time since start = " << total_time << "\n\n"; - - time_file << it << " " << step_time << " " << total_time << " " - << compute_time << " " << refine_time << " " << plot_time << "\n"; - time_file.flush(); - } - - std::cout << "NuFI simulation finished in " << total_time << " seconds.\n"; -} - -NuFISolver::NuFISolver() : order(Parameters::FE_DEGREE), poisson(order) { - std::cout << "Initializing dealii Poisson Solver\n"; - poisson.initialize(); -} +NuFISolver::NuFISolver() = default;