diff --git a/libnufi_lib.a b/libnufi_lib.a index b9e417c..65a365c 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/cells.h b/nufi/cells.h index 19ee0c9..9451c41 100644 --- a/nufi/cells.h +++ b/nufi/cells.h @@ -15,7 +15,7 @@ using namespace dealii; template struct CellInfo { // what needs to be given to evaluator typename DoFHandler::active_cell_iterator cell; - // usefull for locator + // for locator Point lower; Point upper; double h; diff --git a/nufi/parameters.h b/nufi/parameters.h index 2110275..2b0439f 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -31,6 +31,11 @@ constexpr unsigned int FE_DEGREE = 3; constexpr unsigned int CONVERGENCE_ITERATIONS = 5000; constexpr double CONVERGENCE_LIMIT = 1e-8; +// Adaptive refinement options +constexpr unsigned int REFINE_FREQUENCY = 30; +constexpr double REFINEMENT_TOP_FRACTION = 0.8; +constexpr double REFINEMENT_BOTTOM_FRACTION = 0.1; + // Gauge options constexpr double GAUGE_DOMAIN_LEFT = 3.2; constexpr double GAUGE_DOMAIN_RIGHT = 3.8; @@ -42,7 +47,6 @@ constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi) // NUFI options constexpr double DT = 1. / 10.; constexpr unsigned int TMAX = 100; -constexpr unsigned int REFINE_FREQUENCY = 30; // Plotting options constexpr int PLOT_FREQUENCY = 10; diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h index e8abb85..e8548fe 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -1,9 +1,7 @@ #ifndef POISSON_PROBLEM_H #define POISSON_PROBLEM_H -#include #include - #include #include #include @@ -42,6 +40,7 @@ #include #include +#include #include #include #include @@ -49,6 +48,7 @@ #include #include #include +#include #include #include #include @@ -91,6 +91,7 @@ public: return constraints; } const CellLocator &get_locator() const { return cell_locator; } + double get_error_estimate() const { return error_estimate; } std::vector sample_electric_field(double x_min, double x_max, unsigned int Nx); @@ -108,6 +109,7 @@ private: void setup_system(); void assemble_system(); void solve(size_t it); + void estimate_error(); std::function(const std::vector> &)> rhs_function; @@ -124,6 +126,7 @@ private: Vector system_rhs; const bool PRINT_GAUGE_DOF_POSITION = true; + double error_estimate = 0.0; }; //====//====// @@ -336,10 +339,6 @@ template void PoissonProblem::setup_system() { // used for evaluator to avoid running it anytime there is an eval cell_locator.rebuild(dof_handler, triangulation); - - // local_solution_buffer.resize(fe.n_dofs_per_cell()); - // evaluator = std::make_unique>(mapping, fe, - // update_gradients); } template void PoissonProblem::assemble_system() { @@ -419,21 +418,15 @@ template void PoissonProblem::coarse_and_refine_grid(size_t it) { dof_handler, QGauss(fe.degree + 1), std::map *>(), solution, error_per_cell); - GridRefinement::refine_and_coarsen_fixed_number(triangulation, error_per_cell, - 0.3, 0.03); - // START: remove refinment flags from edges of domain to alow safe gauge - // fixing - // for (const auto &cell : triangulation.active_cell_iterators()) { - // const double x = cell->center()[0]; - // if (x >= Parameters::X_DOMAIN_RIGHT - .5) { - // cell->clear_refine_flag(); - // cell->clear_coarsen_flag(); - // } - // } - // END + // GridRefinement::refine_and_coarsen_fixed_number(triangulation, + // error_per_cell, + // 0.3, 0.03); + GridRefinement::refine_and_coarsen_fixed_fraction( + triangulation, error_per_cell, Parameters::REFINEMENT_TOP_FRACTION, + Parameters::REFINEMENT_BOTTOM_FRACTION, + std::numeric_limits::max(), VectorTools::L2_norm); - // triangulation.prepare_coarsening_and_refinement(); triangulation.execute_coarsening_and_refinement(); std::cout << "Refinement Finished..." << "\n"; @@ -450,6 +443,17 @@ template void PoissonProblem::coarse_and_refine_grid(size_t it) { // save_space_vector(Ex, "Ex_after_coarsed", it); } +template void PoissonProblem::estimate_error() { + Vector error_per_cell(triangulation.n_active_cells()); + + KellyErrorEstimator::estimate( + dof_handler, QGauss(fe.degree + 1), + std::map *>(), solution, + error_per_cell); + + error_estimate = error_per_cell.l2_norm(); +} + template void PoissonProblem::solve(size_t it) { std::cout << "Calling PoissonProblem::solve for time-step " << it << "\n"; @@ -499,6 +503,7 @@ void PoissonProblem::solve_step( } assemble_system(); solve(it); + estimate_error(); } // NuFI doesnt use this, kept only for testing PoissonProblem diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index b0044af..c876e28 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -191,27 +191,15 @@ void NuFISolver::run() { using std::abs; using std::max; - // std::unique_ptr rho{ - // reinterpret_cast(std::aligned_alloc(64, sizeof(double) * - // Nx)), - // std::free}; - // - // if (rho == nullptr) - // throw std::bad_alloc{}; - std::vector int_E_squared; int_E_squared.reserve(Nt); std::vector> grid_versions; std::vector> phi_history; - // update_grid_versions(grid_versions, poisson); - // update_solution_history(phi_history, poisson, - // grid_versions.back().grid_version); - std::vector x_eval(Parameters::CALC_NX); - std::ofstream time_file("results/simulation_time.dat"); + std::ofstream time_file(Parameters::PLOT_DIR + "simulation_time.dat"); double total_time = 0; stopwatch total_timer; @@ -224,6 +212,9 @@ void NuFISolver::run() { << "plot_time" << "\n"; + std::ofstream error_file(Parameters::PLOT_DIR + "error_estimate.dat"); + error_file << "it error_estimate\n"; + [[maybe_unused]] const double x_min = Parameters::X_DOMAIN_LEFT; [[maybe_unused]] double dx = Parameters::CALC_DX; @@ -295,6 +286,9 @@ void NuFISolver::run() { 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;