issue source not found, solver correct but breaks after refinement

nufi/grids.h

@@ -2,11 +2,13 @@
 #define GRIDS_H
 
 #include "nufi/cells.h"
+#include <deal.II/base/exceptions.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 <deal.II/numerics/vector_tools.h>
 #include <memory>
 #include <vector>
 
@@ -25,7 +27,7 @@ template <int dim> struct GridStructure {
 
   CellLocator<dim> locator;
 
-  unsigned int grid_version;
+  unsigned int grid_version = 0;
 
   // === // === //
   // Evaluator  //
@@ -35,49 +37,68 @@ template <int dim> struct GridStructure {
                    const std::vector<Point<dim>> &points) const {
 
     std::vector<double> values(points.size());
-#pragma omp parallel
+
-    {
+#pragma omp parallel for
-      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];
-      }
     }
+
     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> grid;
-  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());
+
+  grid.dof_handler = std::make_unique<DoFHandler<dim>>(*grid.triangulation);
+
+  grid.dof_handler->distribute_dofs(*grid.fe);
+
   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;
 }
 

nufi/parameters.h

@@ -35,10 +35,10 @@ 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 = 10;
+constexpr unsigned int REFINE_FREQUENCY = 3;
 
 // Plotting options
-constexpr int PLOT_FREQUENCY = 10;
+constexpr int PLOT_FREQUENCY = 1;
 constexpr size_t PLOT_NX = CALC_NX;
 constexpr double PLOT_DX = LX / PLOT_NX;
 const std::string PLOT_DIR = "results/";

src/nufi_solver.cc

@@ -191,8 +191,8 @@ void NuFISolver::run() {
   std::vector<SolutionSnapshot<1>> phi_history;
 
   update_grid_versions(grid_versions, poisson);
-  update_solution_history(phi_history, poisson,
+  // update_solution_history(phi_history, poisson,
-                          grid_versions.back().grid_version);
+  //                         grid_versions.back().grid_version);
 
   std::vector<double> x_eval(Parameters::CALC_NX);
 
@@ -239,7 +239,21 @@ void NuFISolver::run() {
     // }
     //
     // std::cout << "h ratio = " << max_h / min_h << std::endl;
-
+    // std::vector<double> x = make_x_eval(Parameters::CALC_NX);
+    // auto rho = eval_rho(it, x, grid_versions, phi_history, Parameters::NV);
+    //
+    // double mean = 0;
+    //
+    // for (auto r : rho)
+    //   mean += r;
+    //
+    // mean /= rho.size();
+    //
+    // std::cout << "rho mean = " << mean << "\n";
+    // std::cout << "rho min = " << *std::min_element(rho.begin(), rho.end())
+    //           << "\n";
+    // std::cout << "rho max = " << *std::max_element(rho.begin(), rho.end())
+    //           << "\n";
     // END: diagnostics
 
     double compute_start = timer.elapsed();
@@ -277,24 +291,27 @@ void NuFISolver::run() {
     double step_time = timer_elapsed - time_elapsed_before;
 
     std::cout << "step made in " << step_time << " seconds\n\n";
+
+    //====//====//
+    // Plotting //
+    //====//====//
     if (it % Parameters::PLOT_FREQUENCY == 0) {
       double plot_start = timer.elapsed();
+
       std::cout << "Saving results...   ";
       save_f(*this, it, grid_versions, phi_history, Parameters::PLOT_NX,
              Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
       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());
+      std::vector<double> tmp_Ex(x_eval_Ex.size());
-      tmp_rho = 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);
 
       std::vector<double> E_x(Parameters::PLOT_NX);
       for (size_t i = 0; i < Parameters::PLOT_NX; ++i)
-        E_x[i] = -tmp_rho[i];
+        E_x[i] = -tmp_Ex[i];
       save_space_vector(E_x, "field", it);
 
       double int_val = 0.5 * integral_space_vector_squared(