new paradigm, grids saved when changed old solutions not interpolated to new grids

nufi/cells.h

@@ -6,7 +6,9 @@
 #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 <memory>
 #include <vector>
 
 using namespace dealii;

nufi/fields.h

@@ -1,6 +1,7 @@
 #ifndef FIELDS_H
 #define FIELDS_H
 
+#include "grids.h"
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"
 #include <cmath>
@@ -11,54 +12,6 @@
 
 using namespace dealii;
 
-// inline std::vector<int> Indices_of_points(const std::vector<double> &points,
-// double x_min, double x_max, double dx, int grid_type=0)
-// {
-//     // grid type:
-//     // 0 => uniform
-//     // 1 => non uniform (TODO)
-//
-//     if (dx <= 0.0) {
-//         throw std::invalid_argument("dx must be positive");
-//     }
-//     if (x_max <= x_min) {
-//         throw std::invalid_argument("x_max must be > x_min");
-//     }
-//
-//     std::vector<int> indices;
-//     indices.reserve(points.size());
-//
-//     switch (grid_type) {
-//         case 0:
-//             {
-//             const double L = x_max - x_min;
-//             const int N = std::floor(L/dx);
-//
-//
-//             for (double x : points) //GPT loop, to check
-//             {
-//                 x-= x_min;
-//                 x = x - L * std::floor(x/L);
-//
-//                 int i = static_cast<int>(std::floor(x / dx));
-//
-//                 // safety: handle rare edge case due to floating precision
-//                 if (i == N) i = 0;
-//
-//                 indices.push_back(i);
-//             }
-//             }
-//         case 1:
-//             {
-//             throw std::invalid_argument("Case for non uniform grid is not
-//             completed");
-//             }
-//         default:
-//             throw std::invalid_argument("Invalid grid_type argument");
-//
-//     }
-//     return indices;
-// }
 inline std::vector<double> make_x_eval(size_t Nx) {
   std::vector<double> x_eval(Nx);
   const double dx = Parameters::LX / Nx;
@@ -86,7 +39,7 @@ inline double f0(const double x, const double v,
 
 // wrapper for eval_point() { VectorTools::point_values() }
 inline std::vector<double> eval(std::vector<double> &X,
-                                const PoissonProblem<1> &poisson,
+                                const GridStructure<1> &grid,
                                 const Vector<double> &solution) noexcept {
   size_t x_size = X.size();
   std::vector<double> evals(x_size);
@@ -99,10 +52,10 @@ inline std::vector<double> eval(std::vector<double> &X,
     Points[i][0] = X[i];
   }
 
-  return poisson.eval_vector_grad(solution, Points);
+  return grid.eval_vector_grad(solution, Points);
 }
 
-inline double integral_space_vector(const PoissonProblem<1> &poisson,
+inline double integral_space_vector(const GridStructure<1> &grid,
                                     const Vector<double> &solution,
                                     double dx = Parameters::PLOT_DX,
                                     size_t Nx = Parameters::PLOT_NX) {
@@ -112,13 +65,13 @@ inline double integral_space_vector(const PoissonProblem<1> &poisson,
   for (size_t i = 0; i < Nx; ++i)
     x_eval[i] = xmin + i * dx;
 
-  std::vector<double> tmp = eval(x_eval, poisson, solution);
+  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 PoissonProblem<1> &poisson,
+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) {
@@ -128,7 +81,7 @@ inline double integral_space_vector_squared(const PoissonProblem<1> &poisson,
   for (size_t i = 0; i < Nx; ++i)
     x_eval[i] = xmin + i * dx;
 
-  std::vector<double> tmp = eval(x_eval, poisson, solution);
+  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;

nufi/grids.h

@@ -0,0 +1,116 @@
+#ifndef GRIDS_H
+#define GRIDS_H
+
+#include "nufi/cells.h"
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+#include <memory>
+#include <vector>
+
+using namespace dealii;
+
+template <int dim> class PoissonProblem;
+
+template <int dim> struct GridStructure {
+  //==//==//
+  // Vars //
+  //==//==//
+  std::unique_ptr<Triangulation<dim>> triangulation;
+  std::unique_ptr<DoFHandler<dim>> dof_handler;
+  std::unique_ptr<MappingQ<dim>> mapping;
+  std::unique_ptr<FE_Q<dim>> fe;
+
+  CellLocator<dim> locator;
+
+  unsigned int grid_version;
+
+  // === // === //
+  // Evaluator  //
+  // === // === //
+  std::vector<double>
+  eval_vector_grad(const Vector<double> &solution,
+                   const std::vector<Point<dim>> &points) const {
+
+    std::vector<double> values(points.size());
+#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];
+      }
+    }
+    return values;
+  }
+};
+
+template <int dim>
+GridStructure<dim> make_grid_snapshot(const PoissonProblem<dim> &poisson) {
+  GridStructure<dim> grid;
+  grid.grid_version = 0;
+
+  grid.triangulation = std::make_unique<Triangulation<dim>>();
+  grid.triangulation->copy_triangulation(poisson.get_triangulation());
+
+  grid.mapping = std::make_unique<MappingQ<dim>>(poisson.get_mapping());
+
+  grid.dof_handler = std::make_unique<DoFHandler<dim>>(*grid.triangulation);
+  grid.dof_handler->distribute_dofs(poisson.get_dof_handler().get_fe());
+
+  grid.locator.rebuild(*grid.dof_handler, *grid.triangulation);
+
+  grid.fe = std::make_unique<FE_Q<dim>>(poisson.get_fe());
+
+  return grid;
+}
+
+template <int dim> struct SolutionSnapshot {
+  unsigned int grid_version;
+  Vector<double> solution;
+};
+
+template <int dim>
+inline void update_grid_versions(std::vector<GridStructure<dim>> &grid_versions,
+                                 PoissonProblem<dim> &poisson) {
+  auto grid = make_grid_snapshot(poisson);
+
+  if (!grid_versions.empty())
+    grid.grid_version = grid_versions.back().grid_version + 1;
+
+  grid_versions.push_back(std::move(grid));
+}
+
+template <int dim>
+inline void
+update_solution_history(std::vector<SolutionSnapshot<dim>> &solution_history,
+                        PoissonProblem<dim> &poisson,
+                        unsigned int current_grid_version) {
+
+  SolutionSnapshot<dim> snapshot;
+
+  snapshot.grid_version = current_grid_version;
+
+  Vector<double> solution = poisson.get_solution();
+  snapshot.solution = solution;
+
+  solution_history.push_back(std::move(snapshot));
+}
+
+#endif // !GRIDS_H

nufi/nufi_solver.h

@@ -9,6 +9,7 @@
 #include <vector>
 
 #include "nufi/fields.h" //dont remove
+#include "nufi/grids.h"
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"
 
@@ -19,18 +20,19 @@ public:
   NuFISolver();
 
   void run();
-  std::vector<double> eval_rho(unsigned int n, std::vector<double> &x,
-                               const PoissonProblem<1> &poisson,
-                               const std::vector<Vector<double>> &phi_history,
+  std::vector<double>
+  eval_rho(unsigned int n, std::vector<double> &x,
+           const std::vector<GridStructure<1>> &grid_struct,
+           const std::vector<SolutionSnapshot<1>> &phi_history,
            const unsigned int Nv = Parameters::NV) const;
   std::vector<double>
   eval_ftilda(unsigned int, std::vector<double> &x, double u,
-              const PoissonProblem<1> &poisson,
-              const std::vector<Vector<double>> &phi_history) const;
+              const std::vector<GridStructure<1>> &grid_struct,
+              const std::vector<SolutionSnapshot<1>> &phi_history) const;
   std::vector<double>
   eval_f(unsigned int n, std::vector<double> &x, double u,
-         const PoissonProblem<1> &poisson,
-         const std::vector<Vector<double>> &phi_history) const;
+         const std::vector<GridStructure<1>> &grid_struct,
+         const std::vector<SolutionSnapshot<1>> &phi_history) const;
 
 private:
   unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT);

nufi/poisson_problem.h

@@ -54,7 +54,9 @@
 #include <vector>
 
 #include "nufi/cells.h"
+#include "nufi/grids.h"
 #include "nufi/parameters.h"
+#include "omp.h"
 
 using namespace dealii;
 
@@ -66,8 +68,7 @@ public:
 
   void initialize();
   void solve_step();
-  void coarse_and_refine_grid(size_t it,
-                              std::vector<Vector<double>> &solution_history);
+  void coarse_and_refine_grid(size_t it);
   void run();
 
   unsigned int get_rhs_size();
@@ -78,6 +79,9 @@ public:
   const Vector<double> &get_solution() const { return solution; }
   const MappingQ<dim> &get_mapping() const { return mapping; }
   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 CellLocator<dim> &get_locator() const { return cell_locator; }
 
   std::vector<double> sample_electric_field(double x_min, double x_max,
                                             unsigned int Nx);
@@ -117,8 +121,8 @@ private:
 
   MappingQ<dim> mapping;
 
-  mutable std::vector<double> local_solution_buffer;
-  mutable std::unique_ptr<FEPointEvaluation<dim, dim>> evaluator;
+  // mutable std::vector<double> local_solution_buffer;
+  // mutable std::unique_ptr<FEPointEvaluation<dim, dim>> evaluator;
 };
 
 //====//====//
@@ -211,31 +215,6 @@ PoissonProblem<dim>::sample_electric_potential(double x_min, double x_max,
   return values;
 }
 
-template <int dim>
-std::vector<double> PoissonProblem<dim>::eval_vector_grad(
-    const Vector<double> &solution,
-    const std::vector<Point<dim>> &points) const {
-
-  std::vector<double> values(points.size());
-
-  for (unsigned int p = 0; p < points.size(); ++p) {
-
-    const auto cell_location = cell_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];
-  }
-
-  return values;
-}
-
 template <int dim>
 std::vector<double>
 eval_point_grad(const Mapping<dim> &mapping, const DoFHandler<dim> &dof_handler,
@@ -347,9 +326,9 @@ 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);
+  // local_solution_buffer.resize(fe.n_dofs_per_cell());
+  // evaluator = std::make_unique<FEPointEvaluation<dim, dim>>(mapping, fe,
+  // update_gradients);
 }
 
 // Paul
@@ -402,9 +381,7 @@ template <int dim> void PoissonProblem<dim>::assemble_system() {
   }
 }
 
-template <int dim>
-void PoissonProblem<dim>::coarse_and_refine_grid(
-    size_t it, std::vector<Vector<double>> &solution_history) {
+template <int dim> void PoissonProblem<dim>::coarse_and_refine_grid(size_t it) {
   std::cout << "Refinement Started" << "\n";
   Vector<float> error_per_cell(triangulation.n_active_cells());
 
@@ -417,29 +394,20 @@ void PoissonProblem<dim>::coarse_and_refine_grid(
                                                   0.3, 0.03);
 
   triangulation.prepare_coarsening_and_refinement();
+
   SolutionTransfer<dim, Vector<double>> transfer(dof_handler);
-  transfer.prepare_for_coarsening_and_refinement(solution_history);
+  const Vector<double> refined_solution = solution;
+
+  transfer.prepare_for_coarsening_and_refinement(refined_solution);
 
   triangulation.execute_coarsening_and_refinement();
 
   setup_system();
 
-  std::vector<Vector<double>> new_solution_history(solution_history.size());
-
-  for (auto &vec : new_solution_history)
-    vec.reinit(dof_handler.n_dofs());
-
-  transfer.interpolate(solution_history, new_solution_history);
-
-  solution_history.swap(new_solution_history);
-
-  solution = solution_history.back();
+  transfer.interpolate(refined_solution, solution);
 
   constraints.distribute(solution);
 
-  // cell_locator.rebuild(dof_handler, triangulation); // No need to be called
-  // again because its in setup_system();
-
   std::cout << "Refinement Finished" << "\n";
 
   std::string grid_file_name =

nufi/save_results.h

@@ -1,23 +1,25 @@
 #ifndef SAVE_RESULTS_H
 #define SAVE_RESULTS_H
 
+#include "nufi/grids.h"
 #include "nufi/nufi_solver.h"
 #include "nufi/poisson_problem.h"
 #include <deal.II/lac/vector.h>
 #include <string>
+#include <vector>
 
 void save_f(const NuFISolver &solver, unsigned int n,
-            const PoissonProblem<1> &poisson,
-            const std::vector<Vector<double>> &phi_history, unsigned int Nx_out,
+            std::vector<GridStructure<1>> &grid_struct,
+            std::vector<SolutionSnapshot<1>> &phi_history, unsigned int Nx_out,
             unsigned int Nv_out, const std::string &filename);
 
 void save_rho(const NuFISolver &solver, unsigned int n,
-              const PoissonProblem<1> &poisson,
-              const std::vector<Vector<double>> &phi_history,
+              std::vector<GridStructure<1>> &grid_struct,
+              std::vector<SolutionSnapshot<1>> &phi_history,
               unsigned int Nx_out, const std::string &filename);
 
-void save_Efield(unsigned int n, const PoissonProblem<1> &poisson,
-                 const std::vector<Vector<double>> &phi_history,
+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_space_vector(const std::vector<double> &vals,

src/nufi_solver.cc

@@ -17,6 +17,7 @@
 #include <vector>
 
 #include "nufi/fields.h"
+#include "nufi/grids.h"
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"
 #include "nufi/save_results.h"
@@ -24,10 +25,10 @@
 
 using namespace dealii;
 
-std::vector<double>
-NuFISolver::eval_ftilda(unsigned int n, std::vector<double> &X, double u,
-                        const PoissonProblem<1> &poisson,
-                        const std::vector<Vector<double>> &phi_history) const {
+std::vector<double> NuFISolver::eval_ftilda(
+    unsigned int n, std::vector<double> &X, double u,
+    const std::vector<GridStructure<1>> &grid_struct,
+    const std::vector<SolutionSnapshot<1>> &phi_history) const {
 
   size_t x_size = X.size();
 
@@ -48,7 +49,8 @@ NuFISolver::eval_ftilda(unsigned int n, std::vector<double> &X, double u,
     for (size_t i = 0; i < x_size; ++i)
       X[i] = X[i] - Parameters::DT * U[i];
 
-    tmp = eval(X, poisson, phi_history[n]); // call eval only once
+    tmp = eval(X, grid_struct[phi_history[n].grid_version],
+               phi_history[n].solution); // call eval only once
 
     for (size_t i = 0; i < x_size; ++i) {
       Ex[i] = -tmp[i];
@@ -60,7 +62,8 @@ NuFISolver::eval_ftilda(unsigned int n, std::vector<double> &X, double u,
   for (size_t i = 0; i < x_size; ++i)
     X[i] = X[i] - Parameters::DT * U[i];
 
-  tmp = eval(X, poisson, phi_history[n]); // call eval only once
+  tmp = eval(X, grid_struct[phi_history[n].grid_version],
+             phi_history[n].solution); // call eval only once
 
   for (size_t i = 0; i < x_size; ++i) {
     Ex[i] = -tmp[i];
@@ -74,8 +77,8 @@ NuFISolver::eval_ftilda(unsigned int n, std::vector<double> &X, double u,
 
 std::vector<double>
 NuFISolver::eval_f(unsigned int n, std::vector<double> &X, double u,
-                   const PoissonProblem<1> &poisson,
-                   const std::vector<Vector<double>> &phi_history) const {
+                   const std::vector<GridStructure<1>> &grid_struct,
+                   const std::vector<SolutionSnapshot<1>> &phi_history) const {
 
   size_t x_size = X.size();
 
@@ -92,7 +95,8 @@ NuFISolver::eval_f(unsigned int n, std::vector<double> &X, double u,
 
   std::vector<double> tmp(x_size);
   // Initial half-step.
-  tmp = eval(X, poisson, phi_history[n]); // call eval only once
+  tmp = eval(X, grid_struct[phi_history[n].grid_version],
+             phi_history[n].solution); // call eval only once
   for (size_t i = 0; i < x_size; ++i) {
     Ex[i] = -tmp[i];
     U[i] = U[i] + 0.5 * Parameters::DT * Ex[i];
@@ -102,7 +106,8 @@ NuFISolver::eval_f(unsigned int n, std::vector<double> &X, double u,
     for (size_t i = 0; i < x_size; ++i)
       X[i] = X[i] - Parameters::DT * U[i];
 
-    tmp = eval(X, poisson, phi_history[n]); // call eval only once
+    tmp = eval(X, grid_struct[phi_history[n].grid_version],
+               phi_history[n].solution); // call eval only once
     for (size_t i = 0; i < x_size; ++i) {
       Ex[i] = -tmp[i];
       U[i] = U[i] + Parameters::DT * Ex[i];
@@ -113,7 +118,8 @@ NuFISolver::eval_f(unsigned int n, std::vector<double> &X, double u,
   for (size_t i = 0; i < x_size; ++i)
     X[i] = X[i] - Parameters::DT * U[i];
 
-  tmp = eval(X, poisson, phi_history[n]); // call eval only once
+  tmp = eval(X, grid_struct[phi_history[n].grid_version],
+             phi_history[n].solution); // call eval only once
   for (size_t i = 0; i < x_size; ++i) {
     Ex[i] = -tmp[i];
     U[i] = U[i] + 0.5 * Parameters::DT * Ex[i];
@@ -127,20 +133,21 @@ NuFISolver::eval_f(unsigned int n, std::vector<double> &X, double u,
 
 std::vector<double>
 NuFISolver::eval_rho(unsigned int n, std::vector<double> &X,
-                     const PoissonProblem<1> &poisson,
-                     const std::vector<Vector<double>> &phi_history,
+                     const std::vector<GridStructure<1>> &grid_struct,
+                     const std::vector<SolutionSnapshot<1>> &phi_history,
                      const unsigned int Nv) const {
   size_t x_size = X.size();
+
   const double dv =
       (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
+
   const double v_min = Parameters::V_DOMAIN_LEFT + 0.5 * dv;
 
   std::vector<double> integral(x_size, 0.0);
-
   std::vector<double> tmp_int(x_size);
 
   for (unsigned int i = 0; i < Nv; ++i) {
-    tmp_int = eval_ftilda(n, X, v_min + i * dv, poisson,
+    tmp_int = eval_ftilda(n, X, v_min + i * dv, grid_struct,
                           phi_history); // used eval_ftilda once per i
     for (size_t ii = 0; ii < x_size; ++ii)
       integral[ii] += tmp_int[ii];
@@ -164,16 +171,21 @@ void NuFISolver::run() {
       reinterpret_cast<double *>(std::aligned_alloc(64, sizeof(double) * Nx)),
       std::free};
 
+  if (rho == nullptr)
+    throw std::bad_alloc{};
+
   std::vector<double> int_E_squared;
   int_E_squared.reserve(Nt);
 
-  std::vector<Vector<double>> phi_history;
+  std::vector<GridStructure<1>> grid_versions;
+  std::vector<SolutionSnapshot<1>> phi_history;
+
+  update_grid_versions(grid_versions, poisson);
+  update_solution_history(phi_history, poisson,
+                          grid_versions.back().grid_version);
 
   std::vector<double> x_eval(Parameters::CALC_NX);
 
-  if (rho == nullptr)
-    throw std::bad_alloc{};
-
   std::ofstream time_file("results/simulation_time.dat");
 
   double total_time = 0;
@@ -224,7 +236,7 @@ void NuFISolver::run() {
     // compute rho
     std::vector<double> x_eval = make_x_eval(poisson.get_dof_size());
     std::vector<double> rho_values =
-        eval_rho(it, x_eval, poisson, phi_history, Parameters::NV);
+        eval_rho(it, x_eval, grid_versions, phi_history, Parameters::NV);
 
     Vector<double> rhs(x_eval.size());
     for (unsigned int i = 0; i < rhs.size(); ++i)
@@ -233,18 +245,23 @@ void NuFISolver::run() {
     poisson.set_rhs(rhs);
 
     poisson.solve_step();
-    phi_history.push_back(poisson.get_solution());
 
     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, phi_history);
+
+      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);
 
     double timer_elapsed = timer.elapsed();
     double step_time = timer_elapsed - time_elapsed_before;
@@ -253,30 +270,33 @@ void NuFISolver::run() {
     if (it % Parameters::PLOT_FREQUENCY == 0) {
       double plot_start = timer.elapsed();
       std::cout << "Saving results...   ";
-      save_f(*this, it, poisson, phi_history, Parameters::PLOT_NX,
+      save_f(*this, it, grid_versions, phi_history, Parameters::PLOT_NX,
              Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
-      save_rho(*this, it, poisson, phi_history, Parameters::PLOT_NX,
+      save_rho(*this, it, grid_versions, phi_history, Parameters::PLOT_NX,
                "results/rho_" + std::to_string(it) + ".dat");
       // save_Efield(it, coeffs.get(), 128, "results/field_" +
       // std::to_string(it) + ".dat");
 
       std::vector<double> x_eval_Ex = make_x_eval(Parameters::PLOT_NX);
       std::vector<double> tmp_rho(x_eval_Ex.size());
-      tmp_rho = eval(x_eval_Ex, poisson, phi_history[it]);
+      tmp_rho = eval(x_eval_Ex, grid_versions[phi_history[it].grid_version],
+                     phi_history[it].solution);
 
       std::vector<double> E_x(Parameters::PLOT_NX);
       for (size_t i = 0; i < Parameters::PLOT_NX; ++i)
         E_x[i] = -tmp_rho[i];
       save_space_vector(E_x, "field", it);
 
-      double int_val =
-          0.5 * integral_space_vector_squared(poisson, phi_history[it]);
+      double int_val = 0.5 * integral_space_vector_squared(
+                                 grid_versions[phi_history[it].grid_version],
+                                 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;
     }
+    total_time = timer.elapsed();
 
     time_file << it << " " << step_time << " " << total_time << " "
               << compute_time << " " << refine_time << " " << plot_time << "\n";

src/save_results.cc

@@ -1,18 +1,20 @@
 #include "nufi/save_results.h"
 
 #include "nufi/fields.h"
+#include "nufi/grids.h"
 #include "nufi/nufi_solver.h"
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"
 #include <cstddef>
+#include <deal.II/numerics/solution_transfer.h>
 #include <fstream>
 #include <stdexcept>
 #include <string>
 #include <vector>
 
 void save_f(const NuFISolver &solver, unsigned int n,
-            const PoissonProblem<1> &poisson,
-            const std::vector<Vector<double>> &phi_history, unsigned int Nx_out,
+            std::vector<GridStructure<1>> &grid_struct,
+            std::vector<SolutionSnapshot<1>> &phi_history, unsigned int Nx_out,
             unsigned int Nv_out, const std::string &filename) {
   std::ofstream file(filename);
 
@@ -33,7 +35,7 @@ void save_f(const NuFISolver &solver, unsigned int n,
 
   for (unsigned int j = 0; j < Nv_out; ++j) {
     double v = vmin + (j + 0.5) * dv;
-    val = solver.eval_f(n, x_eval, v, poisson, phi_history);
+    val = solver.eval_f(n, x_eval, v, grid_struct, phi_history);
 
     for (unsigned int i = 0; i < Nx_out; ++i) {
       file << val[i];
@@ -49,8 +51,8 @@ void save_f(const NuFISolver &solver, unsigned int n,
 }
 
 void save_rho(const NuFISolver &solver, unsigned int n,
-              const PoissonProblem<1> &poisson,
-              const std::vector<Vector<double>> &phi_history,
+              std::vector<GridStructure<1>> &grid_struct,
+              std::vector<SolutionSnapshot<1>> &phi_history,
               unsigned int Nx_out, const std::string &filename) {
   std::ofstream file(filename);
 
@@ -61,7 +63,8 @@ void save_rho(const NuFISolver &solver, unsigned int n,
   file << Nx_out << "\n";
   file << xmin << " " << xmax << "\n";
 
-  std::vector<double> tmp = solver.eval_rho(n, x_eval, poisson, phi_history);
+  std::vector<double> tmp =
+      solver.eval_rho(n, x_eval, grid_struct, phi_history);
   for (size_t i = 0; i < Nx_out; ++i) {
     file << tmp[i];
     file << "\n";
@@ -69,9 +72,8 @@ void save_rho(const NuFISolver &solver, unsigned int n,
   file.close();
 }
 
-void save_Efield([[maybe_unused]] unsigned int n,
-                 const PoissonProblem<1> &poisson,
-                 const std::vector<Vector<double>> &phi_history,
+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) {
   std::ofstream file(filename);
 
@@ -85,7 +87,7 @@ void save_Efield([[maybe_unused]] unsigned int n,
   file << Nx_out << "\n";
   file << xmin << " " << xmax << "\n";
 
-  std::vector<double> tmp = eval(x_eval, poisson, phi_history[n]);
+  std::vector<double> tmp = eval(x_eval, grid_struct, phi_history[n].solution);
   for (size_t i = 0; i < Nx_out; ++i) {
     file << -tmp[i];
     file << "\n";