diff --git a/libnufi_lib.a b/libnufi_lib.a index 34b9fa3..8cf9535 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/cells.h b/nufi/cells.h index f30c239..19ee0c9 100644 --- a/nufi/cells.h +++ b/nufi/cells.h @@ -8,7 +8,6 @@ #include #include #include -#include #include using namespace dealii; diff --git a/nufi/fields.h b/nufi/fields.h index 7cb1f83..0d75a8d 100644 --- a/nufi/fields.h +++ b/nufi/fields.h @@ -13,11 +13,19 @@ using namespace dealii; inline std::vector make_x_eval(size_t Nx) { - std::vector x_eval(Nx); - const double dx = Parameters::LX / Nx; - for (size_t i = 0; i < Nx; ++i) - x_eval[i] = Parameters::X_DOMAIN_LEFT + i * dx; - return x_eval; + + std::vector x_eval_E; + + double dx = (Parameters::X_DOMAIN_RIGHT - Parameters::X_DOMAIN_LEFT) / Nx; + + for (unsigned int i = 0; i < Nx; ++i) { + x_eval_E.push_back(Parameters::X_DOMAIN_LEFT + (i + 0.5) * dx); + } + // std::vector x_eval(Nx); + // const double dx = Parameters::LX / Nx; + // for (size_t i = 0; i < Nx; ++i) + // x_eval[i] = Parameters::X_DOMAIN_LEFT + i * dx; + return x_eval_E; } inline void reset_x_eval(std::vector &x_vals) { @@ -41,8 +49,11 @@ inline double f0(const double x, const double v, inline std::vector eval(std::vector &X, const GridStructure<1> &grid, const Vector &solution) noexcept { + + AssertThrow(grid.dof_handler->n_dofs() == solution.size(), + ExcMessage("@ eval(...) grid's number of DoFs doesn't correspond " + "to solution's size")); size_t x_size = X.size(); - std::vector evals(x_size); std::vector> Points(x_size); for (size_t i = 0; i < x_size; ++i) { @@ -97,5 +108,4 @@ Point_vector_to_double_vector(const std::vector> &Points) { return vector; } - #endif diff --git a/nufi/grids.h b/nufi/grids.h index 0c97fc5..be7f960 100644 --- a/nufi/grids.h +++ b/nufi/grids.h @@ -2,13 +2,21 @@ #define GRIDS_H #include "nufi/cells.h" +#include "nufi/parameters.h" + #include #include +#include #include #include #include +#include +#include +#include #include +#include #include + #include #include @@ -24,8 +32,10 @@ template struct GridStructure { std::unique_ptr> dof_handler; std::unique_ptr> mapping; std::unique_ptr> fe; + std::unique_ptr> constraints; + std::unique_ptr sparsity_pattern; - CellLocator locator; + std::unique_ptr> locator; unsigned int grid_version = 0; @@ -38,55 +48,55 @@ template struct GridStructure { std::vector values(points.size()); -#pragma omp parallel for - for (unsigned int p = 0; p < points.size(); ++p) { + // #pragma omp parallel for + // for (unsigned int p = 0; p < points.size(); ++p) { + // + // const Tensor<1, dim> grad_phi = VectorTools::point_gradient( + // *mapping, *dof_handler, solution, points[p]); + // + // values[p] = grad_phi[0]; + // } + // + // return values; + // } - const Tensor<1, dim> grad_phi = VectorTools::point_gradient( - *mapping, *dof_handler, solution, points[p]); +#pragma omp parallel + { + std::vector local_solution_buffer(fe->n_dofs_per_cell()); + FEPointEvaluation evaluator(*mapping, *fe, update_gradients); +#pragma omp for + for (unsigned int p = 0; p < points.size(); ++p) { - values[p] = grad_phi[0]; + const auto cell_location = locator->locate(points[p]); + + cell_location.info->cell->get_dof_values(solution, + local_solution_buffer.begin(), + local_solution_buffer.end()); + + evaluator.reinit( + cell_location.info->cell, + ArrayView>(&cell_location.reference_point, 1)); + + evaluator.evaluate(local_solution_buffer, EvaluationFlags::gradients); + + values[p] = evaluator.get_gradient(0)[0]; + } } - + // AssertThrow(dof_handler->n_dofs() == solution.size(), + // ExcMessage("Solution vector size does not match + // DoFHandler.)"); return values; } - // #pragma omp parallel - // { - // std::vector local_solution_buffer(fe->n_dofs_per_cell()); - // FEPointEvaluation evaluator(*mapping, *fe, - // update_gradients); - // #pragma omp for - // for (unsigned int p = 0; p < points.size(); ++p) { - // - // const auto cell_location = locator.locate(points[p]); - // - // cell_location.info->cell->get_dof_values(solution, - // local_solution_buffer.begin(), - // local_solution_buffer.end()); - // - // evaluator.reinit( - // cell_location.info->cell, - // ArrayView>(&cell_location.reference_point, - // 1)); - // - // evaluator.evaluate(local_solution_buffer, - // EvaluationFlags::gradients); - // - // values[p] = evaluator.get_gradient(0)[0]; - // } - // } - // AssertThrow(dof_handler->n_dofs() == solution.size(), - // ExcMessage("Solution vector size does not match - // DoFHandler.")); - // return values; - // } }; template -GridStructure make_grid_snapshot(const PoissonProblem &poisson) { +GridStructure make_grid_snapshot(PoissonProblem &poisson) { + bool PRINT_GAUGE_DOF_POSITION = true; GridStructure grid; - grid.triangulation = std::make_unique>(); + grid.grid_version = 0; + grid.triangulation = std::make_unique>(); grid.triangulation->copy_triangulation(poisson.get_triangulation()); grid.fe = std::make_unique>(poisson.get_fe()); @@ -96,8 +106,95 @@ GridStructure make_grid_snapshot(const PoissonProblem &poisson) { grid.dof_handler->distribute_dofs(*grid.fe); grid.mapping = std::make_unique>(poisson.get_mapping()); + grid.constraints = + std::make_unique>(poisson.get_constraints()); - grid.locator.rebuild(*grid.dof_handler, *grid.triangulation); + grid.locator = std::make_unique>(); + + // START: transfer poisson.solution to saved grid dofs + SolutionTransfer solution_transfer(*grid.dof_handler); + const Vector coarse_solution = poisson.solution; + solution_transfer.prepare_for_coarsening_and_refinement(coarse_solution); + // START: setup_system(); + grid.constraints->clear(); + + DoFTools::make_hanging_node_constraints(*grid.dof_handler, *grid.constraints); + DoFTools::make_periodicity_constraints(*grid.dof_handler, 0, 1, 0, + *grid.constraints); + + const auto support_points = + DoFTools::map_dofs_to_support_points(*grid.mapping, *grid.dof_handler); + + types::global_dof_index gauge_dof = numbers::invalid_dof_index; + + // Search only inside the protected region + for (const auto &[dof, point] : support_points) { + if (grid.constraints->is_constrained(dof)) + continue; + const double x = point[0]; + if (x <= Parameters::X_DOMAIN_LEFT + .5) { + gauge_dof = dof; + + if (PRINT_GAUGE_DOF_POSITION) + std::cout << " gauge_dof = " << gauge_dof + << " gauge_point = " << point[0] << std::endl; + break; + } + } + + Assert(gauge_dof != numbers::invalid_dof_index, + ExcMessage("No gauge DoF found in protected gauge region.")); + + grid.constraints->add_line(gauge_dof); + grid.constraints->set_inhomogeneity(gauge_dof, 0.0); + + grid.constraints->close(); + + DynamicSparsityPattern dsp(grid.dof_handler->n_dofs()); + DoFTools::make_sparsity_pattern(*grid.dof_handler, dsp, *grid.constraints); + poisson.sparsity_pattern.copy_from(dsp); + poisson.system_matrix.reinit(poisson.sparsity_pattern); + poisson.solution.reinit(grid.dof_handler->n_dofs()); + poisson.system_rhs.reinit(grid.dof_handler->n_dofs()); + + grid.locator->rebuild(*grid.dof_handler, *grid.triangulation); + // END + + solution_transfer.interpolate(coarse_solution, poisson.solution); + // END + + // START: diagnostics + AssertThrow( + poisson.get_constraints().n_constraints() == + grid.constraints->n_constraints(), + ExcMessage( + "PoissonProblem constraints doesn't match Snapshot constraints")); + for (auto c1 = poisson.get_dof_handler().begin_active(), + c2 = grid.dof_handler->begin_active(); + c1 != poisson.get_dof_handler().end(); ++c1, ++c2) { + std::vector d1(c1->get_fe().dofs_per_cell); + std::vector d2(c2->get_fe().dofs_per_cell); + + c1->get_dof_indices(d1); + c2->get_dof_indices(d2); + + AssertThrow(d1 == d2, ExcInternalError()); + } + // + // auto support_points = + // DoFTools::map_dofs_to_support_points(*grid.mapping, *grid.dof_handler); + // + // std::cout << "COPY\n"; + // + // for (const auto &[dof, point] : support_points) { + // std::cout << dof << " : " << point[0] << "\n"; + // } + // for (auto cell : grid.dof_handler->active_cell_iterators()) { + // std::cout << "@ GridStructure: " << cell->id() << " " << + // cell->center()[0] + // << '\n'; + // } + // END return grid; } @@ -116,6 +213,10 @@ inline void update_grid_versions(std::vector> &grid_versions, grid.grid_version = grid_versions.back().grid_version + 1; grid_versions.push_back(std::move(grid)); + + std::cout << "@ update_grid_history: " + << "grid version " << grid.grid_version << " size " + << poisson.get_solution().size() << "\n"; } template @@ -128,9 +229,15 @@ update_solution_history(std::vector> &solution_history, snapshot.grid_version = current_grid_version; + // where I might need to change something to pass the correct solution or in + // the correct form Vector solution = poisson.get_solution(); snapshot.solution = solution; + std::cout << "@ update_solution_history: " + << "grid version " << current_grid_version << " size " + << poisson.get_solution().size() << "\n"; + solution_history.push_back(std::move(snapshot)); } diff --git a/nufi/parameters.h b/nufi/parameters.h index 56ee566..675ec72 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -23,14 +23,14 @@ constexpr unsigned int NV = 128; constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV; // deal.ii options -constexpr unsigned int GLOBAL_REFINEMENT = 6; +constexpr unsigned int GLOBAL_REFINEMENT = 8; constexpr unsigned int FE_DEGREE = 3; constexpr unsigned int CONVERGENCE_ITERATIONS = 5000; constexpr double CONVERGENCE_LIMIT = 1e-8; -// gauge fix options -constexpr double GAUGE_X_MIN = 2.5; -constexpr double GAUGE_X_MAX = 3.5; +// Gauge options +constexpr double GAUGE_DOMAIN_LEFT = 3.2; +constexpr double GAUGE_DOMAIN_RIGHT = 3.8; constexpr double EPS = 0.01; constexpr double WAVE_NR = 0.5; @@ -39,7 +39,7 @@ 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 = 1; +constexpr unsigned int REFINE_FREQUENCY = 3; // Plotting options constexpr int PLOT_FREQUENCY = 1; diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h index bf8ba37..3f3dbc9 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -1,6 +1,7 @@ #ifndef POISSON_PROBLEM_H #define POISSON_PROBLEM_H +#include #include #include @@ -69,8 +70,10 @@ public: PoissonProblem(unsigned int degree); void initialize(); - void solve_step(); + void solve_step(size_t it, std::vector> &grid_versions, + bool refining = false); void coarse_and_refine_grid(size_t it); + void setup_constraints(AffineConstraints &constraints); void run(); unsigned int get_rhs_size(); @@ -84,6 +87,9 @@ public: const DoFHandler &get_dof_handler() const { return dof_handler; } const Triangulation &get_triangulation() const { return triangulation; } const FE_Q &get_fe() const { return fe; } + const AffineConstraints &get_constraints() const { + return constraints; + } const CellLocator &get_locator() const { return cell_locator; } std::vector sample_electric_field(double x_min, double x_max, @@ -101,12 +107,16 @@ public: DoFHandler dof_handler; CellLocator cell_locator; + Vector solution; // phi + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector system_rhs; + private: void create_mesh(); void setup_system(); - void setup_gauge(); void assemble_system(); - void solve(); + void solve(size_t it); // Triangulation triangulation; FE_Q fe; @@ -114,11 +124,11 @@ private: AffineConstraints constraints; - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; + // SparsityPattern sparsity_pattern; + // SparseMatrix system_matrix; - Vector solution; // phi - Vector system_rhs; + // Vector solution; // phi + // Vector system_rhs; std::function(const std::vector> &)> rhs_function; @@ -154,8 +164,9 @@ void PoissonProblem::set_rhs_function( template PoissonProblem::PoissonProblem(unsigned int degree) - : triangulation(Triangulation::limit_level_difference_at_vertices), - dof_handler(triangulation), fe(degree), mapping(degree) {} + : // triangulation(Triangulation::limit_level_difference_at_vertices), + triangulation(), dof_handler(triangulation), fe(degree), mapping(degree) { +} template std::vector @@ -291,7 +302,13 @@ template void PoissonProblem::create_mesh() { triangulation.refine_global(Parameters::GLOBAL_REFINEMENT); } -template void PoissonProblem::setup_gauge() { +template +void PoissonProblem::setup_constraints( + AffineConstraints &constraints) { + + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + DoFTools::make_periodicity_constraints(dof_handler, 0, 1, 0, constraints); + const auto support_points = DoFTools::map_dofs_to_support_points(mapping, dof_handler); @@ -302,8 +319,7 @@ template void PoissonProblem::setup_gauge() { if (constraints.is_constrained(dof)) continue; const double x = point[0]; - if (x >= Parameters::X_DOMAIN_RIGHT - .5 || - (x <= Parameters::X_DOMAIN_LEFT + .5 && x != 0)) { + if (x <= Parameters::X_DOMAIN_LEFT + .5) { gauge_dof = dof; if (PRINT_GAUGE_DOF_POSITION) @@ -326,17 +342,14 @@ template void PoissonProblem::setup_system() { constraints.clear(); - DoFTools::make_hanging_node_constraints(dof_handler, constraints); - DoFTools::make_periodicity_constraints(dof_handler, 0, 1, 0, constraints); - - setup_gauge(); + setup_constraints(constraints); constraints.close(); + // DSP DynamicSparsityPattern dsp(dof_handler.n_dofs()); DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); sparsity_pattern.copy_from(dsp); - system_matrix.reinit(sparsity_pattern); solution.reinit(dof_handler.n_dofs()); @@ -418,7 +431,7 @@ template void PoissonProblem::assemble_system() { } template void PoissonProblem::coarse_and_refine_grid(size_t it) { - std::cout << "Refinement Started" << "\n"; + std::cout << "Refinement Started..." << "\n"; Vector error_per_cell(triangulation.n_active_cells()); KellyErrorEstimator::estimate( @@ -432,30 +445,17 @@ template void PoissonProblem::coarse_and_refine_grid(size_t it) { // fixing // for (const auto &cell : triangulation.active_cell_iterators()) { // const double x = cell->center()[0]; - // if (x >= Parameters::X_DOMAIN_RIGHT - .5 || - // x <= Parameters::X_DOMAIN_LEFT + .5) { + // if (x >= Parameters::X_DOMAIN_RIGHT - .5) { // cell->clear_refine_flag(); // cell->clear_coarsen_flag(); // } // } // END - triangulation.prepare_coarsening_and_refinement(); - - SolutionTransfer> transfer(dof_handler); - const Vector refined_solution = solution; - - transfer.prepare_for_coarsening_and_refinement(refined_solution); - + // triangulation.prepare_coarsening_and_refinement(); triangulation.execute_coarsening_and_refinement(); - setup_system(); - - transfer.interpolate(refined_solution, solution); - - constraints.distribute(solution); - - std::cout << "Refinement Finished" << "\n"; + std::cout << "Refinement Finished..." << "\n"; std::string grid_file_name = Parameters::PLOT_DIR + "grid_" + std::to_string(it); @@ -469,20 +469,36 @@ template void PoissonProblem::coarse_and_refine_grid(size_t it) { // save_space_vector(Ex, "Ex_after_coarsed", it); } -template void PoissonProblem::solve() { +template void PoissonProblem::solve(size_t it) { SolverControl solver_control(Parameters::CONVERGENCE_ITERATIONS, Parameters::CONVERGENCE_LIMIT * system_rhs.l2_norm()); SolverCG> solver(solver_control); - // PreconditionSSOR> preconditioner; - // preconditioner.initialize(system_matrix, 1.2); - - // solver.solve(system_matrix, solution, system_rhs, preconditioner); - solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity()); constraints.distribute(solution); + + std::ofstream out("results/phi_after_solve_" + std::to_string(it) + ".dat"); + std::vector> data; + + const auto support = + DoFTools::map_dofs_to_support_points(mapping, dof_handler); + + for (const auto &[dof, p] : support) { + data.emplace_back(p[0], solution[dof]); + } + + std::sort(data.begin(), data.end()); + + for (const auto &[x, value] : data) { + out << x << " " << value << "\n"; + } + + std::vector E_x = + sample_electric_field(Parameters::X_DOMAIN_LEFT, + Parameters::X_DOMAIN_RIGHT, Parameters::PLOT_NX); + save_space_vector(E_x, "E_x_after_solve", it); } template void PoissonProblem::initialize() { @@ -490,9 +506,19 @@ template void PoissonProblem::initialize() { setup_system(); // distribute DoFs and matrices } -template void PoissonProblem::solve_step() { +template +void PoissonProblem::solve_step( + size_t it, std::vector> &grid_versions, bool refining) { + if (refining) { + coarse_and_refine_grid(it); + setup_system(); + } + assemble_system(); - solve(); + solve(it); + + if (refining) + update_grid_versions(grid_versions, *this); } // NuFI doesnt use this, kept only for testing PoissonProblem diff --git a/nufi/save_results.h b/nufi/save_results.h index 73d9ce0..227d351 100644 --- a/nufi/save_results.h +++ b/nufi/save_results.h @@ -23,6 +23,11 @@ void save_Efield(unsigned int n, GridStructure<1> &grid_struct, std::vector> &phi_history, unsigned int Nx_out, const std::string &filename); +void save_Efield_new(unsigned int it, + std::vector> &grid_versions, + std::vector> &phi_history, + unsigned int Nx_out = Parameters::PLOT_NX); + void save_space_vector(const std::vector &vals, const std::string &filename, size_t it); diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index e4f8e67..3c57079 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -58,6 +58,11 @@ std::vector NuFISolver::eval_ftilda( ", expected by grid_struct = " + std::to_string(grid_struct[phi_history[n].grid_version] .dof_handler->n_dofs()))); + AssertThrow( + grid_struct[phi_history[n].grid_version].grid_version == + phi_history[n].grid_version, + ExcMessage( + "grid.grid_version not equal to phi_history[n].grid_version")); tmp = eval(X, grid_struct[phi_history[n].grid_version], phi_history[n].solution); // call eval only once @@ -200,7 +205,7 @@ void NuFISolver::run() { std::vector> grid_versions; std::vector> phi_history; - update_grid_versions(grid_versions, poisson); + // update_grid_versions(grid_versions, poisson); // update_solution_history(phi_history, poisson, // grid_versions.back().grid_version); @@ -278,21 +283,13 @@ void NuFISolver::run() { return eval_rho(it, x, grid_versions, phi_history, Parameters::NV); }); - poisson.solve_step(); - - 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); - - update_grid_versions(grid_versions, poisson); - - refine_time = timer.elapsed() - refine_start; - std::cout << "Refinement step done in " - << std::to_string(std::round(std::floor(refine_time))) << "[s]" - << "\n"; + if (it % Parameters::REFINE_FREQUENCY == 0) { + // if (it == 0) { + poisson.solve_step(it, grid_versions, true); + compute_time = timer.elapsed() - compute_start; + } 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); @@ -314,15 +311,16 @@ void NuFISolver::run() { save_rho(*this, it, grid_versions, phi_history, Parameters::PLOT_NX, "results/rho_" + std::to_string(it) + ".dat"); - std::vector x_eval_Ex = make_x_eval(Parameters::PLOT_NX); - std::vector tmp_Ex(x_eval_Ex.size()); - tmp_Ex = eval(x_eval_Ex, grid_versions[phi_history[it].grid_version], - phi_history[it].solution); - - std::vector E_x(Parameters::PLOT_NX); - for (size_t i = 0; i < Parameters::PLOT_NX; ++i) - E_x[i] = -tmp_Ex[i]; - save_space_vector(E_x, "field", it); + // std::vector x_eval_Ex = make_x_eval(Parameters::PLOT_NX); + // std::vector tmp_Ex(x_eval_Ex.size()); + // tmp_Ex = eval(x_eval_Ex, grid_versions[phi_history[it].grid_version], + // phi_history[it].solution); + // + // std::vector E_x(Parameters::PLOT_NX); + // for (size_t i = 0; i < Parameters::PLOT_NX; ++i) + // E_x[i] = -tmp_Ex[i]; + // save_space_vector(E_x, "field", it); + save_Efield_new(it, grid_versions, phi_history); double int_val = 0.5 * integral_space_vector_squared( grid_versions[phi_history[it].grid_version], diff --git a/src/save_results.cc b/src/save_results.cc index bd787d3..431d0cb 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -72,6 +72,35 @@ void save_rho(const NuFISolver &solver, unsigned int n, file.close(); } +void save_Efield_new(unsigned int it, + std::vector> &grid_versions, + std::vector> &phi_history, + unsigned int Nx_out) { + std::vector x_eval_E = make_x_eval(Nx_out); + + auto grad_phi = eval(x_eval_E, grid_versions[phi_history[it].grid_version], + phi_history[it].solution); + + std::vector> ordered_E; + + for (size_t i = 0; i < x_eval_E.size(); ++i) { + double E = -grad_phi[i]; + ordered_E.push_back({x_eval_E[i], E}); + } + + std::sort(ordered_E.begin(), ordered_E.end()); + + std::ofstream E_file("results/E_" + std::to_string(it) + ".dat"); + + for (auto [x, E] : ordered_E) + E_file << x << " " << E << "\n"; + + std::cout << "Saving iteration " << it << " using grid version " + << phi_history[it].grid_version << "\n"; + + E_file.close(); +} + void save_Efield([[maybe_unused]] unsigned int n, GridStructure<1> &grid_struct, std::vector> &phi_history, unsigned int Nx_out, const std::string &filename) {