diff --git a/libnufi_lib.a b/libnufi_lib.a index dc5e371..34b9fa3 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/parameters.h b/nufi/parameters.h index dc35759..56ee566 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -28,6 +28,10 @@ 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; + constexpr double EPS = 0.01; constexpr double WAVE_NR = 0.5; constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi) diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h index 8cced11..bf8ba37 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -104,8 +104,7 @@ public: private: void create_mesh(); void setup_system(); - void setup_system_in_refinement_before_interpolating(); - void setup_system_in_refinement_after_interpolating(); + void setup_gauge(); void assemble_system(); void solve(); @@ -127,6 +126,8 @@ private: MappingQ mapping; + const bool PRINT_GAUGE_DOF_POSITION = true; + // mutable std::vector local_solution_buffer; // mutable std::unique_ptr> evaluator; }; @@ -290,6 +291,35 @@ template void PoissonProblem::create_mesh() { triangulation.refine_global(Parameters::GLOBAL_REFINEMENT); } +template void PoissonProblem::setup_gauge() { + const auto support_points = + DoFTools::map_dofs_to_support_points(mapping, 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 (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)) { + 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.")); + + constraints.add_line(gauge_dof); + constraints.set_inhomogeneity(gauge_dof, 0.0); +} + template void PoissonProblem::setup_system() { dof_handler.distribute_dofs(fe); @@ -297,26 +327,9 @@ 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); - // Gauge fix for periodic Poisson: - // remove the constant nullspace by pinning one unconstrained DoF. - // (by Paul Wilhelm) - types::global_dof_index gauge_dof = numbers::invalid_dof_index; - - for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i) { - if (!constraints.is_constrained(i)) { - gauge_dof = i; - break; - } - } - - Assert(gauge_dof != numbers::invalid_dof_index, - ExcMessage("No unconstrained DoF found for gauge fixing.")); - - constraints.add_line(gauge_dof); - constraints.set_inhomogeneity(gauge_dof, 0.0); + setup_gauge(); constraints.close(); @@ -337,55 +350,6 @@ template void PoissonProblem::setup_system() { // update_gradients); } -template -void PoissonProblem::setup_system_in_refinement_before_interpolating() { - dof_handler.distribute_dofs(fe); - - constraints.clear(); - - DoFTools::make_hanging_node_constraints(dof_handler, constraints); - - DoFTools::make_periodicity_constraints(dof_handler, 0, 1, 0, constraints); - - // Do NOT close the constraints yet. - // The gauge will be added after interpolation. - - solution.reinit(dof_handler.n_dofs()); - system_rhs.reinit(dof_handler.n_dofs()); - - cell_locator.rebuild(dof_handler, triangulation); -} - -template -void PoissonProblem::setup_system_in_refinement_after_interpolating() { - - // Add the gauge constraint. - types::global_dof_index gauge_dof = numbers::invalid_dof_index; - - for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i) { - if (!constraints.is_constrained(i)) { - gauge_dof = i; - break; - } - } - - Assert(gauge_dof != numbers::invalid_dof_index, - ExcMessage("No unconstrained DoF found for gauge fixing.")); - - constraints.add_line(gauge_dof); - constraints.set_inhomogeneity(gauge_dof, 0.0); - - constraints.close(); - - DynamicSparsityPattern dsp(dof_handler.n_dofs()); - - DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); - - sparsity_pattern.copy_from(dsp); - - system_matrix.reinit(sparsity_pattern); -} - template void PoissonProblem::assemble_system() { Assert(system_matrix.m() == dof_handler.n_dofs(), @@ -464,6 +428,18 @@ template void PoissonProblem::coarse_and_refine_grid(size_t it) { 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 || + // x <= Parameters::X_DOMAIN_LEFT + .5) { + // cell->clear_refine_flag(); + // cell->clear_coarsen_flag(); + // } + // } + // END + triangulation.prepare_coarsening_and_refinement(); SolutionTransfer> transfer(dof_handler); @@ -474,12 +450,9 @@ template void PoissonProblem::coarse_and_refine_grid(size_t it) { triangulation.execute_coarsening_and_refinement(); setup_system(); - // setup_system_in_refinement_before_interpolating(); transfer.interpolate(refined_solution, solution); - // setup_system_in_refinement_after_interpolating(); - constraints.distribute(solution); std::cout << "Refinement Finished" << "\n"; diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index e99acab..e4f8e67 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -329,11 +329,12 @@ void NuFISolver::run() { phi_history[it].solution); 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; + std::cout << "Results saved in " << plot_start << "[s]" << "\n"; } total_time = 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";