diff --git a/libnufi_lib.a b/libnufi_lib.a index 708e8e1..11289ac 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/parameters.h b/nufi/parameters.h index c379718..4fe02a2 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -3,6 +3,7 @@ #include #include +#include namespace Parameters { constexpr unsigned int DIMENSION = 1; @@ -24,7 +25,7 @@ constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV; // deal.ii options constexpr unsigned int GLOBAL_REFINEMENT = 6; constexpr unsigned int FE_DEGREE = 2; -constexpr unsigned int CONVERGENCE_ITERATIONS = 5000; +constexpr unsigned int CONVERGENCE_ITERATIONS = 10000; constexpr double CONVERGENCE_LIMIT = 1e-8; constexpr double EPS = 0.01; @@ -32,13 +33,15 @@ constexpr double WAVE_NR = 0.5; constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi) // NUFI options -constexpr double DT = 1. / 8.; +constexpr double DT = 1. / 10.; constexpr unsigned int TMAX = 100; +constexpr unsigned int REFINE_FREQUENCY = 5; // Plotting options -constexpr int PLOT_FREQUENCY = 4; +constexpr int PLOT_FREQUENCY = 5; constexpr size_t PLOT_NX = CALC_NX; constexpr double PLOT_DX = LX / PLOT_NX; +const std::string PLOT_DIR = "results/"; } // namespace Parameters #endif diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h index 475cd62..2ca9355 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -66,7 +66,8 @@ public: void initialize(); void solve_step(); - void coarse_and_refine_grid(); + void coarse_and_refine_grid(size_t it, + std::vector> &solution_history); void run(); void set_rhs_function(std::function &)> f); @@ -84,7 +85,10 @@ public: eval_vector_grad(const Vector &solution, const std::vector> &points) const; - void save_grid_to_file(const std::string &filename) const; + void save_grid_to_file(std::string &filename) const; + + Triangulation triangulation; + DoFHandler dof_handler; private: void create_mesh(); @@ -92,9 +96,9 @@ private: void assemble_system(); void solve(); - Triangulation triangulation; + // Triangulation triangulation; FE_Q fe; - DoFHandler dof_handler; + // DoFHandler dof_handler; AffineConstraints constraints; @@ -250,10 +254,18 @@ double eval_point_value(const Mapping &mapping, } template -void PoissonProblem::save_grid_to_file(const std::string &filename) const { - std::ofstream out(filename); +void PoissonProblem::save_grid_to_file(std::string &filename) const { GridOut grid_out; - grid_out.write_svg(triangulation, out); + + if (dim >= 2) { + filename += ".svg"; + std::ofstream out(filename); + grid_out.write_svg(triangulation, out); + } else if (dim == 1) { + filename += ".gnuplot"; + std::ofstream out(filename); + grid_out.write_gnuplot(triangulation, out); + } std::cout << "Grid written to " << filename << "\n"; } @@ -371,12 +383,8 @@ template void PoissonProblem::assemble_system() { template void PoissonProblem::coarse_and_refine_grid( - std::vector> &solution_history) { - // Add refinement and coasring algorithm here - - //======//======// - // CHECK step-15. - //======//======// + size_t it, std::vector> &solution_history) { + std::cout << "Refinement Started" << "\n"; Vector error_per_cell(triangulation.n_active_cells()); KellyErrorEstimator::estimate( @@ -384,15 +392,10 @@ void PoissonProblem::coarse_and_refine_grid( std::map *>(), solution, error_per_cell); - GridRefinement::refine_and_coarsen_fixed_number( - triangulation, error_per_cell, 0.3, - 0.03); // These numbers are to be reviewd and changed + GridRefinement::refine_and_coarsen_fixed_number(triangulation, error_per_cell, + 0.3, 0.03); triangulation.prepare_coarsening_and_refinement(); - - // old solutions transfer - const size_t N_sols = solution_history.size(); - SolutionTransfer> transfer(dof_handler); transfer.prepare_for_coarsening_and_refinement(solution_history); @@ -400,14 +403,27 @@ void PoissonProblem::coarse_and_refine_grid( setup_system(); - transfer.interpolate(solution_history); + std::vector> new_solution_history(solution_history.size()); + + for (auto &vec : new_solution_history) + vec.reinit(dof_handler.n_dofs()); + + transfer.interpolate(solution_history, new_solution_history); + + solution_history.swap(new_solution_history); solution = solution_history.back(); constraints.distribute(solution); - // rebuild cells with the grid cell_locator.rebuild(dof_handler, triangulation); + + std::cout << "Refinement Finished" << "\n"; + + std::string grid_file_name = + Parameters::PLOT_DIR + "grid_" + std::to_string(it); + + save_grid_to_file(grid_file_name); } template void PoissonProblem::solve() { diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index 75637eb..ee95a31 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -173,7 +173,14 @@ void NuFISolver::run() { std::ofstream time_file("results/simulation_time.dat"); double total_time = 0; - time_file << "it step_time total_time" << "\n"; + + time_file << "it " + << "step_time " + << "total_time " + << "compute_time " + << "refine_time " + << "plot_time" + << "\n"; const double x_min = Parameters::X_DOMAIN_LEFT; double dx = Parameters::CALC_DX; @@ -182,11 +189,30 @@ void NuFISolver::run() { 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 << ")" << std::endl; + // START: diagnostics + std::cout << "cells = " << poisson.triangulation.n_active_cells() + << " dofs = " << poisson.dof_handler.n_dofs() << std::endl; + double min_h = 1e100; + double max_h = 0; + + for (auto cell : poisson.triangulation.active_cell_iterators()) { + min_h = std::min(min_h, cell->diameter()); + max_h = std::max(max_h, cell->diameter()); + } + + std::cout << "h ratio = " << max_h / min_h << std::endl; + + // END: diagnostics + + double compute_start = timer.elapsed(); // compute rho std::vector x_eval = make_x_eval(Nx); std::vector tmp_rho = @@ -207,19 +233,20 @@ void NuFISolver::run() { poisson.solve_step(); phi_history.push_back(poisson.get_solution()); - // std::vector sampled_potential = - // poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of - // FE + compute_time = timer.elapsed() - compute_start; + + if (it % Parameters::REFINE_FREQUENCY == 0 && it != 0) { + double refine_start = timer.elapsed(); + poisson.coarse_and_refine_grid(it, phi_history); + refine_time = timer.elapsed() - refine_start; + } double timer_elapsed = timer.elapsed(); double step_time = timer_elapsed - time_elapsed_before; - total_time += timer_elapsed; - - time_file << it << " " << step_time << " " << total_time << "\n"; - time_file.flush(); 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... "; save_f(*this, it, poisson, phi_history, Parameters::PLOT_NX, Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat"); @@ -241,7 +268,13 @@ void NuFISolver::run() { int_E_squared.push_back(int_val); save_space_vector(int_E_squared, "electricint", it); std::cout << "Time since start = " << total_time << "\n\n"; + + plot_time = timer.elapsed() - plot_start; } + + 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";