more trials on refinement bug, tried with grid update rebuiltind setup constraints

nufi/cells.h

@@ -8,7 +8,6 @@
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/fe/mapping_q.h>
 #include <deal.II/grid/tria.h>
-#include <memory>
 #include <vector>
 
 using namespace dealii;

nufi/fields.h

@@ -13,11 +13,19 @@
 using namespace dealii;
 
 inline std::vector<double> make_x_eval(size_t Nx) {
-  std::vector<double> x_eval(Nx);
+
-  const double dx = Parameters::LX / Nx;
+  std::vector<double> x_eval_E;
-  for (size_t i = 0; i < Nx; ++i)
+
-    x_eval[i] = Parameters::X_DOMAIN_LEFT + i * dx;
+  double dx = (Parameters::X_DOMAIN_RIGHT - Parameters::X_DOMAIN_LEFT) / Nx;
-  return x_eval;
+
+  for (unsigned int i = 0; i < Nx; ++i) {
+    x_eval_E.push_back(Parameters::X_DOMAIN_LEFT + (i + 0.5) * dx);
+  }
+  // std::vector<double> 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<double> &x_vals) {
@@ -41,8 +49,11 @@ inline double f0(const double x, const double v,
 inline std::vector<double> eval(std::vector<double> &X,
                                 const GridStructure<1> &grid,
                                 const Vector<double> &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<double> evals(x_size);
   std::vector<Point<1>> Points(x_size);
 
   for (size_t i = 0; i < x_size; ++i) {
@@ -97,5 +108,4 @@ Point_vector_to_double_vector(const std::vector<Point<1>> &Points) {
 
   return vector;
 }
-
 #endif

nufi/grids.h

@@ -2,13 +2,21 @@
 #define GRIDS_H
 
 #include "nufi/cells.h"
+#include "nufi/parameters.h"
+
 #include <deal.II/base/exceptions.h>
 #include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/mapping_q.h>
 #include <deal.II/grid/tria.h>
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
 #include <deal.II/matrix_free/fe_point_evaluation.h>
+#include <deal.II/numerics/solution_transfer.h>
 #include <deal.II/numerics/vector_tools.h>
+
 #include <memory>
 #include <vector>
 
@@ -24,8 +32,10 @@ template <int dim> struct GridStructure {
   std::unique_ptr<DoFHandler<dim>> dof_handler;
   std::unique_ptr<MappingQ<dim>> mapping;
   std::unique_ptr<FE_Q<dim>> fe;
+  std::unique_ptr<AffineConstraints<double>> constraints;
+  std::unique_ptr<SparsityPattern> sparsity_pattern;
 
-  CellLocator<dim> locator;
+  std::unique_ptr<CellLocator<dim>> locator;
 
   unsigned int grid_version = 0;
 
@@ -38,55 +48,55 @@ template <int dim> struct GridStructure {
 
     std::vector<double> values(points.size());
 
-#pragma omp parallel for
+    // #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;
-  }
-  // #pragma omp parallel
-  //     {
-  //       std::vector<double> local_solution_buffer(fe->n_dofs_per_cell());
-  //       FEPointEvaluation<dim, dim> 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]);
+    //       const Tensor<1, dim> grad_phi = VectorTools::point_gradient(
+    //           *mapping, *dof_handler, solution, points[p]);
     //
-  //         cell_location.info->cell->get_dof_values(solution,
+    //       values[p] = grad_phi[0];
-  //                                                  local_solution_buffer.begin(),
-  //                                                  local_solution_buffer.end());
-  //
-  //         evaluator.reinit(
-  //             cell_location.info->cell,
-  //             ArrayView<const Point<dim>>(&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<double> local_solution_buffer(fe->n_dofs_per_cell());
+      FEPointEvaluation<dim, dim> 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<const Point<dim>>(&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 <int dim>
-GridStructure<dim> make_grid_snapshot(const PoissonProblem<dim> &poisson) {
+GridStructure<dim> make_grid_snapshot(PoissonProblem<dim> &poisson) {
+  bool PRINT_GAUGE_DOF_POSITION = true;
   GridStructure<dim> grid;
 
-  grid.triangulation = std::make_unique<Triangulation<dim>>();
+  grid.grid_version = 0;
 
+  grid.triangulation = std::make_unique<Triangulation<dim>>();
   grid.triangulation->copy_triangulation(poisson.get_triangulation());
 
   grid.fe = std::make_unique<FE_Q<dim>>(poisson.get_fe());
@@ -96,8 +106,95 @@ GridStructure<dim> make_grid_snapshot(const PoissonProblem<dim> &poisson) {
   grid.dof_handler->distribute_dofs(*grid.fe);
 
   grid.mapping = std::make_unique<MappingQ<dim>>(poisson.get_mapping());
+  grid.constraints =
+      std::make_unique<AffineConstraints<double>>(poisson.get_constraints());
 
-  grid.locator.rebuild(*grid.dof_handler, *grid.triangulation);
+  grid.locator = std::make_unique<CellLocator<dim>>();
+
+  // START: transfer poisson.solution to saved grid dofs
+  SolutionTransfer<dim> solution_transfer(*grid.dof_handler);
+  const Vector<double> 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<types::global_dof_index> d1(c1->get_fe().dofs_per_cell);
+    std::vector<types::global_dof_index> 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<GridStructure<dim>> &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 <int dim>
@@ -128,9 +229,15 @@ update_solution_history(std::vector<SolutionSnapshot<dim>> &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<double> 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));
 }
 

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
+// Gauge options
-constexpr double GAUGE_X_MIN = 2.5;
+constexpr double GAUGE_DOMAIN_LEFT = 3.2;
-constexpr double GAUGE_X_MAX = 3.5;
+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;

nufi/poisson_problem.h

@@ -1,6 +1,7 @@
 #ifndef POISSON_PROBLEM_H
 #define POISSON_PROBLEM_H
 
+#include <cstddef>
 #include <deal.II/base/function.h>
 
 #include <deal.II/base/index_set.h>
@@ -69,8 +70,10 @@ public:
   PoissonProblem(unsigned int degree);
 
   void initialize();
-  void solve_step();
+  void solve_step(size_t it, std::vector<GridStructure<1>> &grid_versions,
+                  bool refining = false);
   void coarse_and_refine_grid(size_t it);
+  void setup_constraints(AffineConstraints<double> &constraints);
   void run();
 
   unsigned int get_rhs_size();
@@ -84,6 +87,9 @@ public:
   const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
   const Triangulation<dim> &get_triangulation() const { return triangulation; }
   const FE_Q<dim> &get_fe() const { return fe; }
+  const AffineConstraints<double> &get_constraints() const {
+    return constraints;
+  }
   const CellLocator<dim> &get_locator() const { return cell_locator; }
 
   std::vector<double> sample_electric_field(double x_min, double x_max,
@@ -101,12 +107,16 @@ public:
   DoFHandler<dim> dof_handler;
   CellLocator<dim> cell_locator;
 
+  Vector<double> solution; // phi
+  SparsityPattern sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+  Vector<double> system_rhs;
+
 private:
   void create_mesh();
   void setup_system();
-  void setup_gauge();
   void assemble_system();
-  void solve();
+  void solve(size_t it);
 
   // Triangulation<dim> triangulation;
   FE_Q<dim> fe;
@@ -114,11 +124,11 @@ private:
 
   AffineConstraints<double> constraints;
 
-  SparsityPattern sparsity_pattern;
+  // SparsityPattern sparsity_pattern;
-  SparseMatrix<double> system_matrix;
+  // SparseMatrix<double> system_matrix;
 
-  Vector<double> solution; // phi
+  // Vector<double> solution; // phi
-  Vector<double> system_rhs;
+  // Vector<double> system_rhs;
 
   std::function<std::vector<double>(const std::vector<Point<dim>> &)>
       rhs_function;
@@ -154,8 +164,9 @@ void PoissonProblem<dim>::set_rhs_function(
 
 template <int dim>
 PoissonProblem<dim>::PoissonProblem(unsigned int degree)
-    : triangulation(Triangulation<dim>::limit_level_difference_at_vertices),
+    : // triangulation(Triangulation<dim>::limit_level_difference_at_vertices),
-      dof_handler(triangulation), fe(degree), mapping(degree) {}
+      triangulation(), dof_handler(triangulation), fe(degree), mapping(degree) {
+}
 
 template <int dim>
 std::vector<double>
@@ -291,7 +302,13 @@ template <int dim> void PoissonProblem<dim>::create_mesh() {
   triangulation.refine_global(Parameters::GLOBAL_REFINEMENT);
 }
 
-template <int dim> void PoissonProblem<dim>::setup_gauge() {
+template <int dim>
+void PoissonProblem<dim>::setup_constraints(
+    AffineConstraints<double> &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 <int dim> void PoissonProblem<dim>::setup_gauge() {
     if (constraints.is_constrained(dof))
       continue;
     const double x = point[0];
-    if (x >= Parameters::X_DOMAIN_RIGHT - .5 ||
+    if (x <= Parameters::X_DOMAIN_LEFT + .5) {
-        (x <= Parameters::X_DOMAIN_LEFT + .5 && x != 0)) {
       gauge_dof = dof;
 
       if (PRINT_GAUGE_DOF_POSITION)
@@ -326,17 +342,14 @@ template <int dim> void PoissonProblem<dim>::setup_system() {
 
   constraints.clear();
 
-  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+  setup_constraints(constraints);
-  DoFTools::make_periodicity_constraints(dof_handler, 0, 1, 0, constraints);
-
-  setup_gauge();
 
   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 <int dim> void PoissonProblem<dim>::assemble_system() {
 }
 
 template <int dim> void PoissonProblem<dim>::coarse_and_refine_grid(size_t it) {
-  std::cout << "Refinement Started" << "\n";
+  std::cout << "Refinement Started..." << "\n";
 
   Vector<float> error_per_cell(triangulation.n_active_cells());
   KellyErrorEstimator<dim>::estimate(
@@ -432,30 +445,17 @@ template <int dim> void PoissonProblem<dim>::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 ||
+  //   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();
+  // triangulation.prepare_coarsening_and_refinement();
-
-  SolutionTransfer<dim, Vector<double>> transfer(dof_handler);
-  const Vector<double> refined_solution = solution;
-
-  transfer.prepare_for_coarsening_and_refinement(refined_solution);
-
   triangulation.execute_coarsening_and_refinement();
 
-  setup_system();
+  std::cout << "Refinement Finished..." << "\n";
-
-  transfer.interpolate(refined_solution, solution);
-
-  constraints.distribute(solution);
-
-  std::cout << "Refinement Finished" << "\n";
 
   std::string grid_file_name =
       Parameters::PLOT_DIR + "grid_" + std::to_string(it);
@@ -469,20 +469,36 @@ template <int dim> void PoissonProblem<dim>::coarse_and_refine_grid(size_t it) {
   // save_space_vector(Ex, "Ex_after_coarsed", it);
 }
 
-template <int dim> void PoissonProblem<dim>::solve() {
+template <int dim> void PoissonProblem<dim>::solve(size_t it) {
 
   SolverControl solver_control(Parameters::CONVERGENCE_ITERATIONS,
                                Parameters::CONVERGENCE_LIMIT *
                                    system_rhs.l2_norm());
   SolverCG<Vector<double>> solver(solver_control);
 
-  // PreconditionSSOR<SparseMatrix<double>> 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<std::pair<double, double>> 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<double> 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 <int dim> void PoissonProblem<dim>::initialize() {
@@ -490,9 +506,19 @@ template <int dim> void PoissonProblem<dim>::initialize() {
   setup_system(); // distribute DoFs and matrices
 }
 
-template <int dim> void PoissonProblem<dim>::solve_step() {
+template <int dim>
+void PoissonProblem<dim>::solve_step(
+    size_t it, std::vector<GridStructure<1>> &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

nufi/save_results.h

@@ -23,6 +23,11 @@ void save_Efield(unsigned int n, GridStructure<1> &grid_struct,
                  std::vector<SolutionSnapshot<1>> &phi_history,
                  unsigned int Nx_out, const std::string &filename);
 
+void save_Efield_new(unsigned int it,
+                     std::vector<GridStructure<1>> &grid_versions,
+                     std::vector<SolutionSnapshot<1>> &phi_history,
+                     unsigned int Nx_out = Parameters::PLOT_NX);
+
 void save_space_vector(const std::vector<double> &vals,
                        const std::string &filename, size_t it);
 

src/nufi_solver.cc

@@ -58,6 +58,11 @@ std::vector<double> 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<GridStructure<1>> grid_versions;
   std::vector<SolutionSnapshot<1>> 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();
+    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;
-
-    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";
     }
     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<double> x_eval_Ex = make_x_eval(Parameters::PLOT_NX);
+      // std::vector<double> x_eval_Ex = make_x_eval(Parameters::PLOT_NX);
-      std::vector<double> tmp_Ex(x_eval_Ex.size());
+      // std::vector<double> tmp_Ex(x_eval_Ex.size());
-      tmp_Ex = eval(x_eval_Ex, grid_versions[phi_history[it].grid_version],
+      // tmp_Ex = eval(x_eval_Ex, grid_versions[phi_history[it].grid_version],
-                    phi_history[it].solution);
+      //               phi_history[it].solution);
-
+      //
-      std::vector<double> E_x(Parameters::PLOT_NX);
+      // std::vector<double> E_x(Parameters::PLOT_NX);
-      for (size_t i = 0; i < Parameters::PLOT_NX; ++i)
+      // for (size_t i = 0; i < Parameters::PLOT_NX; ++i)
-        E_x[i] = -tmp_Ex[i];
+      //   E_x[i] = -tmp_Ex[i];
-      save_space_vector(E_x, "field", it);
+      // 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],

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<GridStructure<1>> &grid_versions,
+                     std::vector<SolutionSnapshot<1>> &phi_history,
+                     unsigned int Nx_out) {
+  std::vector<double> 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<std::pair<double, double>> 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<SolutionSnapshot<1>> &phi_history,
                  unsigned int Nx_out, const std::string &filename) {