Merge branch 'refinement_tree_locator' into 1x2v

# Conflicts:
#	Makefile
#	README.md
#	cmake_install.cmake
#	libnufi_lib.a
#	nufi/fields.h
#	nufi/parameters.h
#	nufi/save_results.h
#	src/nufi_solver.cc
#	src/save_results.cc

nufi/cells.h

@@ -1,117 +1,146 @@
 #ifndef CELLS_H
 #define CELLS_H
 
-#include <algorithm>
-#include <boost/geometry/geometries/concepts/point_concept.hpp>
+#include <array>
+#include <cmath>
 #include <deal.II/base/geometry_info.h>
 #include <deal.II/base/point.h>
 #include <deal.II/dofs/dof_handler.h>
-#include <deal.II/fe/mapping_q.h>
 #include <deal.II/grid/tria.h>
 #include <vector>
 
+#include "nufi/parameters.h"
+
 using namespace dealii;
 
-template <int dim> struct CellInfo {
-  // what needs to be given to evaluator
-  typename DoFHandler<dim>::active_cell_iterator cell;
-  // usefull for locator
-  Point<dim> lower;
-  Point<dim> upper;
-  double h;
-};
-
 template <int dim> struct CellLocation {
-  const CellInfo<dim> *info;
+  typename DoFHandler<dim>::active_cell_iterator cell;
   Point<dim> reference_point;
 };
 
 template <int dim> class CellLocator {
 public:
-  using CellIterator = typename DoFHandler<dim>::active_cell_iterator;
-
   void rebuild(const DoFHandler<dim> &dof_handler,
                const Triangulation<dim> &triangulation);
+
   CellLocation<dim> locate(const Point<dim> &p) const;
 
-  const std::vector<Point<dim>> &get_cell_centers() const;
-
 private:
-  std::vector<CellInfo<dim>> cells;
-  std::vector<Point<dim>> cell_centers;
+  const DoFHandler<dim> *dof_handler_ptr = nullptr;
+  Point<dim> lower;
+  Point<dim> upper;
+
+  // Cached base level = Parameters::GLOBAL_REFINEMENT.
+  unsigned int base_level = 0;
+  unsigned int base_n_per_axis = 1; // 2^base_level
+
+  std::vector<typename Triangulation<dim>::cell_iterator> base_cells;
 };
 
 template <int dim>
 void CellLocator<dim>::rebuild(const DoFHandler<dim> &dof_handler,
                                const Triangulation<dim> &triangulation) {
+  dof_handler_ptr = &dof_handler;
 
-  cells.clear();
-  cells.reserve(triangulation.n_active_cells());
+  AssertThrow(triangulation.n_cells(0) == 1,
+              ExcMessage("CellLocator assumes exactly one coarse/root cell."));
 
-  for (const auto &cell : dof_handler.active_cell_iterators()) {
-    CellInfo<dim> info;
+  typename Triangulation<dim>::cell_iterator root = triangulation.begin(0);
+  lower = root->vertex(0);
+  upper = root->vertex(GeometryInfo<dim>::vertices_per_cell - 1);
 
-    info.cell = cell;
-    info.lower = cell->vertex(0);
-    info.upper = cell->vertex(GeometryInfo<dim>::vertices_per_cell - 1);
+  base_level = Parameters::GLOBAL_REFINEMENT;
+  base_n_per_axis = 1u << base_level; // = 2^base_level
 
-    info.h = info.upper[0] - info.lower[0];
+  const unsigned int n_base_cells =
+      1u << (dim * base_level); // = 2^(dim*base_level)
+  base_cells.assign(n_base_cells, typename Triangulation<dim>::cell_iterator());
 
-    cells.push_back(info);
-  }
+  std::vector<typename Triangulation<dim>::cell_iterator> stack;
+  std::vector<unsigned int> index_stack;
+  std::vector<unsigned int> depth_stack;
+  stack.push_back(root);
+  index_stack.push_back(0);
+  depth_stack.push_back(0);
 
-  std::sort(cells.begin(), cells.end(),
-            [](const CellInfo<dim> &a, const CellInfo<dim> &b) {
-              return a.lower[0] < b.lower[0];
-            });
+  while (!stack.empty()) {
+    auto cell = stack.back();
+    unsigned int idx = index_stack.back();
+    unsigned int depth = depth_stack.back();
+    stack.pop_back();
+    index_stack.pop_back();
+    depth_stack.pop_back();
 
-  cell_centers.clear();
-  cell_centers.reserve(cells.size());
+    if (depth == base_level) {
+      base_cells[idx] = cell;
+      continue;
+    }
 
-  for (const auto &cell : cells) {
-    Point<dim> center;
-    for (unsigned int d = 0; d < dim; ++d)
-      center[d] = 0.5 * (cell.lower[d] + cell.upper[d]);
+    AssertThrow(cell->has_children(),
+                ExcMessage("CellLocator: mesh is not uniformly refined to "
+                           "Parameters::GLOBAL_REFINEMENT; base-level cache "
+                           "cannot be built. Did you coarsen below the "
+                           "global refinement level?"));
 
-    cell_centers.push_back(center);
+    const unsigned int n_children = GeometryInfo<dim>::max_children_per_cell;
+    for (unsigned int c = 0; c < n_children; ++c) {
+      stack.push_back(cell->child(c));
+      index_stack.push_back(idx * n_children + c);
+      depth_stack.push_back(depth + 1);
+    }
   }
 }
 
 template <int dim>
 CellLocation<dim> CellLocator<dim>::locate(const Point<dim> &p) const {
-  static_assert(dim == 1,
-                "Current CellLocator implementation only supports 1D.");
-
-  AssertThrow(!cells.empty(),
+  AssertThrow(dof_handler_ptr != nullptr,
               ExcMessage("CellLocator::rebuild() has not been called."));
 
-  const double x = p[0];
+  Point<dim> p_wrapped;
+  for (unsigned int d = 0; d < dim; ++d) {
+    const double L = upper[d] - lower[d];
+    double x = p[d] - lower[d];
+    x = x - L * std::floor(x / L);
+    p_wrapped[d] = lower[d] + x;
+  }
 
-  auto it = std::upper_bound(cells.begin(), cells.end(), x,
-                             [](double value, const CellInfo<dim> &cell) {
-                               return value < cell.lower[0];
-                             }); // returns cell to the right of cell with x
+  Point<dim> xi;
+  for (unsigned int d = 0; d < dim; ++d) {
+    xi[d] = (p_wrapped[d] - lower[d]) / (upper[d] - lower[d]);
+    xi[d] = std::min(std::max(xi[d], 0.0), 1.0); // of cell at base level
+  }
 
-  if (it == cells.begin())
-    it = cells.begin();
-  else
-    --it;
+  unsigned int idx = 0;
+  Point<dim> xi_local = xi;
+  for (unsigned int l = 0; l < base_level; ++l) {
+    const unsigned int child_index =
+        GeometryInfo<dim>::child_cell_from_point(xi_local);
+    xi_local =
+        GeometryInfo<dim>::cell_to_child_coordinates(xi_local, child_index);
+    idx = idx * GeometryInfo<dim>::max_children_per_cell + child_index;
+  }
 
-  // Safety check: make sure the point is really inside this cell
-  AssertThrow(x >= it->lower[0] - 1e-12 && x <= it->upper[0] + 1e-12,
-              ExcMessage("CellLocator failed to find containing cell."));
+  typename Triangulation<dim>::cell_iterator cell = base_cells[idx];
+  xi = xi_local;
+
+  // Step 4: continue descending only through ADAPTIVE refinement beyond
+  // the base level -- this loop now only runs `depth - base_level` times
+  // instead of `depth` times.
+  while (cell->has_children()) {
+    const unsigned int child_index =
+        GeometryInfo<dim>::child_cell_from_point(xi);
+    xi = GeometryInfo<dim>::cell_to_child_coordinates(xi, child_index);
+    cell = cell->child(child_index);
+  }
+
+  typename DoFHandler<dim>::active_cell_iterator dof_cell(
+      &cell->get_triangulation(), cell->level(), cell->index(),
+      dof_handler_ptr);
 
   CellLocation<dim> location;
-
-  location.info = &(*it);
-  location.reference_point[0] = (p[0] - it->lower[0]) / it->h;
-
+  location.cell = dof_cell;
+  location.reference_point = xi;
   return location;
 }
 
-template <int dim>
-const std::vector<Point<dim>> &CellLocator<dim>::get_cell_centers() const {
-  return cell_centers;
-}
-
 #endif // !CELLS_H

nufi/fields.h

@@ -1,9 +1,9 @@
-#ifndef FIELDS_H
-#define FIELDS_H
+#ifndef NUFI_FIELDS_H_
+#define NUFI_FIELDS_H_
 
 #include "nufi/grids.h"
 #include "nufi/parameters.h"
-#include "nufi/poisson_problem.h"
+
 #include <cmath>
 #include <cstddef>
 #include <deal.II/base/function.h>
@@ -13,25 +13,76 @@
 using namespace dealii;
 
 inline std::vector<double> make_x_eval(size_t Nx) {
-
   std::vector<double> x_eval_E;
 
   double dx = (Parameters::X_DOMAIN_RIGHT - Parameters::X_DOMAIN_LEFT) / Nx;
 
-  for (unsigned int i = 0; i < Nx; ++i) {
+  for (unsigned int i = 0; i < Nx; ++i)
     x_eval_E.push_back(Parameters::X_DOMAIN_LEFT + (i + 0.5) * dx);
-  }
+
   return x_eval_E;
 }
 
-inline double f0(const double x, const double v1, const double v2,
-                 const double eps = Parameters::EPS,
-                 const double k = Parameters::WAVE_NR) {
-  const double prefactor =
-      Parameters::F0_FACTOR * (1.0 + eps * std::cos(k * x));
-  const double gaussian = v1 * v1 * std::exp(-0.5 * (v1 * v1 + v2 * v2));
+inline void reset_x_eval(std::vector<double> &x_vals) {
+  const size_t Nx = x_vals.size();
+  const double dx = Parameters::LX / Nx;
+  for (size_t i = 0; i < Nx; ++i)
+    x_vals[i] = Parameters::X_DOMAIN_LEFT + i * dx;
+};
 
-  return prefactor * gaussian;
+inline double f0(const double x, const double v,
+                 const size_t f0_type = Parameters::f0_TYPE) {
+  const double eps = Parameters::EPS;
+  const double k = Parameters::WAVE_NR;
+
+  const double factor = Parameters::F0_FACTOR;
+
+  auto maxwell = [factor](double u, double v = 0, double v_th = 1) {
+    return 1 / v_th * factor *
+           std::exp(-0.5 * (u - v) * (u - v) / (v_th * v_th));
+  };
+
+  double prefactor;
+  double computed_max;
+  double result;
+
+  switch (f0_type) {
+  case 0: // two-stream
+  {
+    computed_max = maxwell(v);
+    prefactor = (1.0 + eps * std::cos(k * x)) * v * v;
+    result = prefactor * computed_max;
+    break;
+  }
+  case 1: // Landau-damping
+  {
+    computed_max = maxwell(v);
+    prefactor = (1.0 + eps * std::cos(k * x));
+    result = prefactor * computed_max;
+    break;
+  }
+  case 2: // Maxwellian
+  {
+    result = maxwell(v);
+    break;
+  }
+  case 3: // Bump-on tail
+  {
+    const double beam_density = 0.05;
+    const double beam_v_th = 0.2;
+    const double beam_v = 3;
+    const double alpha = beam_density / (1 - beam_density);
+
+    const double beam_max = maxwell(v, beam_v, beam_v_th);
+    computed_max = maxwell(v);
+    result = (1 - alpha) * computed_max + alpha * beam_max;
+    break;
+  }
+  default:
+    throw std::invalid_argument("Invalid f0_type");
+  }
+
+  return result;
 }
 
 // wrapper for eval_point() { VectorTools::point_values() }
@@ -42,50 +93,41 @@ inline std::vector<double> eval(std::vector<double> &X,
   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<Point<1>> Points(x_size);
 
-  for (size_t i = 0; i < x_size; ++i) {
-    X[i] = X[i] - Parameters::X_DOMAIN_LEFT;
-    X[i] = X[i] - Parameters::LX * std::floor(X[i] * Parameters::LX_INV);
+  const size_t x_size = X.size();
+  std::vector<Point<1>> points(x_size);
 
-    Points[i][0] = X[i];
-  }
+  for (size_t i = 0; i < x_size; ++i)
+    points[i][0] = X[i];
 
-  return grid.eval_vector_grad(solution, Points);
+  return grid.eval_vector_grad(solution, points);
 }
 
 inline double integral_space_vector(const GridStructure<1> &grid,
                                     const Vector<double> &solution,
-                                    double dx = Parameters::PLOT_DX,
                                     size_t Nx = Parameters::PLOT_NX) {
   double integral = 0.0;
-  double xmin = Parameters::X_DOMAIN_LEFT;
-  std::vector<double> x_eval(Nx);
-  for (size_t i = 0; i < Nx; ++i)
-    x_eval[i] = xmin + i * dx;
+  std::vector<double> x_eval = make_x_eval(Nx);
+  const double dx = Parameters::LX / Nx;
 
   std::vector<double> tmp = eval(x_eval, grid, solution);
   for (size_t i = 0; i < Nx; ++i)
     integral += tmp[i];
   return integral * dx;
-};
+}
 
 inline double integral_space_vector_squared(const GridStructure<1> &grid,
                                             const Vector<double> &solution,
-                                            double dx = Parameters::PLOT_DX,
                                             size_t Nx = Parameters::PLOT_NX) {
   double integral = 0.0;
-  double xmin = Parameters::X_DOMAIN_LEFT;
-  std::vector<double> x_eval(Nx);
-  for (size_t i = 0; i < Nx; ++i)
-    x_eval[i] = xmin + i * dx;
+  std::vector<double> x_eval = make_x_eval(Nx);
+  const double dx = Parameters::LX / Nx;
 
   std::vector<double> tmp = eval(x_eval, grid, solution);
   for (size_t i = 0; i < Nx; ++i)
     integral += tmp[i] * tmp[i];
   return integral * dx;
-};
+}
 
 inline std::vector<double>
 Point_vector_to_double_vector(const std::vector<Point<1>> &Points) {
@@ -97,4 +139,4 @@ Point_vector_to_double_vector(const std::vector<Point<1>> &Points) {
 
   return vector;
 }
-#endif
+#endif // NUFI_FIELDS_H_

nufi/grids.h

@@ -72,12 +72,12 @@ template <int dim> struct GridStructure {
 
         const auto cell_location = locator->locate(points[p]);
 
-        cell_location.info->cell->get_dof_values(solution,
-                                                 local_solution_buffer.begin(),
-                                                 local_solution_buffer.end());
+        cell_location.cell->get_dof_values(solution,
+                                           local_solution_buffer.begin(),
+                                           local_solution_buffer.end());
 
         evaluator.reinit(
-            cell_location.info->cell,
+            cell_location.cell,
             ArrayView<const Point<dim>>(&cell_location.reference_point, 1));
 
         evaluator.evaluate(local_solution_buffer, EvaluationFlags::gradients);
@@ -127,7 +127,7 @@ GridStructure<dim> make_grid_snapshot(PoissonProblem<dim> &poisson) {
 
       if (PRINT_GAUGE_DOF_POSITION)
         std::cout << " gauge_dof = " << gauge_dof
-                  << " gauge_point = " << point[0] << std::endl;
+                  << " gauge_point = " << point[0] << "\n";
       break;
     }
   }

nufi/nufi_solver.h

@@ -8,7 +8,7 @@
 #include <deal.II/numerics/vector_tools.h>
 #include <vector>
 
-#include "nufi/fields.h" //dont remove
+#include "nufi/fields.h" // dont remove
 #include "nufi/grids.h"
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"

nufi/parameters.h

@@ -2,6 +2,7 @@
 #define PARAMETERS_H
 
 #include <cmath>
+#include <cstddef>
 #include <cstdlib>
 #include <string>
 
@@ -16,23 +17,30 @@ constexpr double LX_INV = 1 / LX;
 constexpr size_t CALC_NX = 256;
 constexpr double CALC_DX = LX / CALC_NX;
 
-constexpr double V_DOMAIN_LEFT_1 = -10.;
-constexpr double V_DOMAIN_RIGHT_1 = 10.;
+constexpr double V_DOMAIN_LEFT = -10.;
+constexpr double V_DOMAIN_RIGHT = 10.;
 
-constexpr unsigned int NV_1 = 128;
-constexpr double DV_1 = std::abs(V_DOMAIN_RIGHT_1 - V_DOMAIN_LEFT_1) / NV_1;
+constexpr unsigned int NV = 128;
+constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV;
 
-constexpr double V_DOMAIN_LEFT_2 = -10.;
-constexpr double V_DOMAIN_RIGHT_2 = 10.;
+// f0_TYPE:
+// 0 -> twos-stream
+// 1 -> landau-damping
+// 2 -> maxwellian
+// 3 -> bump-on-tail
+constexpr size_t f0_TYPE = 0;
 
-constexpr unsigned int NV_2 = 128;
-constexpr double DV_2 = std::abs(V_DOMAIN_RIGHT_2 - V_DOMAIN_LEFT_2) / NV_2;
 // deal.ii options
-constexpr unsigned int GLOBAL_REFINEMENT = 6;
+constexpr unsigned int GLOBAL_REFINEMENT = 7;
 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;
@@ -44,13 +52,11 @@ 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;
 constexpr size_t PLOT_NX = CALC_NX;
 constexpr double PLOT_DX = LX / PLOT_NX;
-constexpr double PLOT_FIXED_V2 = 0.;
 const std::string PLOT_DIR = "results/";
 } // namespace Parameters
 

nufi/poisson_problem.h

@@ -1,9 +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>
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/mpi_remote_point_evaluation.h>
@@ -42,6 +40,7 @@
 #include <deal.II/numerics/solution_transfer.h>
 #include <deal.II/numerics/vector_tools.h>
 
+#include <cstddef>
 #include <deal.II/numerics/vector_tools_evaluate.h>
 #include <deal.II/numerics/vector_tools_interpolate.h>
 #include <deal.II/numerics/vector_tools_point_gradient.h>
@@ -49,6 +48,7 @@
 #include <fstream>
 #include <functional>
 #include <iostream>
+#include <limits>
 #include <string>
 #include <utility>
 #include <vector>
@@ -56,6 +56,7 @@
 #include "nufi/cells.h"
 #include "nufi/grids.h"
 #include "nufi/parameters.h"
+#include "nufi/stopwatch.h"
 #include "omp.h"
 
 void save_space_vector(const std::vector<double> &vals,
@@ -70,8 +71,8 @@ public:
   PoissonProblem(unsigned int degree);
 
   void initialize();
-  void solve_step(size_t it, std::vector<GridStructure<1>> &grid_versions,
-                  bool refining = false);
+  double 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();
@@ -91,6 +92,7 @@ public:
     return constraints;
   }
   const CellLocator<dim> &get_locator() const { return cell_locator; }
+  double get_error_estimate() const { return error_estimate; }
 
   std::vector<double> sample_electric_field(double x_min, double x_max,
                                             unsigned int Nx);
@@ -108,6 +110,7 @@ private:
   void setup_system();
   void assemble_system();
   void solve(size_t it);
+  void estimate_error();
 
   std::function<std::vector<double>(const std::vector<Point<dim>> &)>
       rhs_function;
@@ -124,6 +127,7 @@ private:
   Vector<double> system_rhs;
 
   const bool PRINT_GAUGE_DOF_POSITION = true;
+  double error_estimate = 0.0;
 };
 
 //====//====//
@@ -336,10 +340,6 @@ template <int dim> void PoissonProblem<dim>::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<FEPointEvaluation<dim, dim>>(mapping, fe,
-  // update_gradients);
 }
 
 template <int dim> void PoissonProblem<dim>::assemble_system() {
@@ -375,7 +375,9 @@ template <int dim> void PoissonProblem<dim>::assemble_system() {
   Assert(rhs_function,
          ExcMessage("Poisson RHS function has not been initialized."));
 
+  std::cout << "Start of full rho eval..." << "\n";
   std::vector<double> all_rho = rhs_function(all_q_points);
+  std::cout << "End of full rho eval..." << "\n";
 
   Assert(all_rho.size() == all_q_points.size(),
          ExcMessage("rhs_function returned wrong size"));
@@ -417,21 +419,20 @@ template <int dim> void PoissonProblem<dim>::coarse_and_refine_grid(size_t it) {
       dof_handler, QGauss<dim - 1>(fe.degree + 1),
       std::map<types::boundary_id, const Function<dim> *>(), 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<unsigned int>::max(), VectorTools::L2_norm);
+
+  // Avoid coarsing below GLOBAL_REFINEMENT level for CellLocator
+  for (const auto &cell : triangulation.active_cell_iterators())
+    if (cell->level() <= static_cast<int>(Parameters::GLOBAL_REFINEMENT))
+      cell->clear_coarsen_flag();
 
-  // triangulation.prepare_coarsening_and_refinement();
   triangulation.execute_coarsening_and_refinement();
 
   std::cout << "Refinement Finished..." << "\n";
@@ -439,13 +440,17 @@ template <int dim> void PoissonProblem<dim>::coarse_and_refine_grid(size_t it) {
   std::string grid_file_name =
       Parameters::PLOT_DIR + "grid_" + std::to_string(it);
   save_grid_to_file(grid_file_name);
+}
 
-  // std::vector<double> Ex =
-  // sample_electric_potential(Parameters::X_DOMAIN_LEFT,
-  //                                                    Parameters::X_DOMAIN_RIGHT,
-  //                                                    Parameters::PLOT_NX);
-  //
-  // save_space_vector(Ex, "Ex_after_coarsed", it);
+template <int dim> void PoissonProblem<dim>::estimate_error() {
+  Vector<float> error_per_cell(triangulation.n_active_cells());
+
+  KellyErrorEstimator<dim>::estimate(
+      dof_handler, QGauss<dim - 1>(fe.degree + 1),
+      std::map<types::boundary_id, const Function<dim> *>(), solution,
+      error_per_cell);
+
+  error_estimate = error_per_cell.l2_norm();
 }
 
 template <int dim> void PoissonProblem<dim>::solve(size_t it) {
@@ -488,15 +493,21 @@ template <int dim> void PoissonProblem<dim>::initialize() {
 }
 
 template <int dim>
-void PoissonProblem<dim>::solve_step(
+double PoissonProblem<dim>::solve_step(
     size_t it, std::vector<GridStructure<1>> &grid_versions, bool refining) {
+  double refining_time = 0.0;
   if (refining) {
+    stopwatch<double> refining_timer;
     coarse_and_refine_grid(it);
     setup_system();
     update_grid_versions(grid_versions, *this);
+    refining_time = refining_timer.elapsed();
   }
   assemble_system();
   solve(it);
+  estimate_error();
+
+  return refining_time;
 }
 
 // NuFI doesnt use this, kept only for testing PoissonProblem

nufi/save_results.h

@@ -8,21 +8,29 @@
 
 class NuFISolver;
 
-void save_f(const NuFISolver &solver, unsigned int n, const double v2_0,
-            std::vector<GridStructure<1>> &grid_struct,
-            std::vector<SolutionSnapshot<1>> &phi_history, unsigned int Nx_out,
-            unsigned int Nv_out, const std::string &filename);
+struct DiagnosticsSnapshot {
+  unsigned int Nx = 0;
+  unsigned int Nv = 0;
+  std::vector<double> x_eval;
+  std::vector<double> v_eval;
+  std::vector<double> f;
+  std::vector<double> rho;
+  std::vector<double> E;
+};
 
-void save_rho(const NuFISolver &solver, unsigned int n,
-              std::vector<GridStructure<1>> &grid_struct,
-              std::vector<SolutionSnapshot<1>> &phi_history,
-              unsigned int Nx_out, const std::string &filename);
+DiagnosticsSnapshot
+compute_diagnostics(const NuFISolver &solver, unsigned int n,
+                    std::vector<GridStructure<1>> &grid_struct,
+                    std::vector<SolutionSnapshot<1>> &phi_history,
+                    unsigned int Nx_out, unsigned int Nv_out);
 
-void save_Efield(unsigned int it, std::vector<GridStructure<1>> &grid_versions,
-                 std::vector<SolutionSnapshot<1>> &phi_history,
-                 unsigned int Nx_out = Parameters::PLOT_NX);
+void save_f(const DiagnosticsSnapshot &snap, const std::string &filepath);
+void save_rho(const DiagnosticsSnapshot &snap, const std::string &filepath);
+void save_Efield(const DiagnosticsSnapshot &snap, const std::string &filepath);
+double compute_int_E_squared(const DiagnosticsSnapshot &snap);
 
-void save_space_vector(const std::vector<double> &vals,
-                       const std::string &filename, size_t it);
+void save_time_series(const std::vector<double> &t,
+                      const std::vector<double> &values,
+                      const std::string &filepath);
 
 #endif