bug identified to be in refinement function, maybe in interpolation step...

nufi/parameters.h

@@ -35,7 +35,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 = 3;
+constexpr unsigned int REFINE_FREQUENCY = 1;
 
 // Plotting options
 constexpr int PLOT_FREQUENCY = 1;

nufi/poisson_problem.h

@@ -56,8 +56,12 @@
 #include "nufi/cells.h"
 #include "nufi/grids.h"
 #include "nufi/parameters.h"
+#include "nufi/save_results.h"
 #include "omp.h"
 
+void save_space_vector(const std::vector<double> &vals,
+                       const std::string &filename, size_t it);
+
 using namespace dealii;
 
 // =-=-=-=-= Poisson Solver =-=-=-=-=
@@ -432,6 +436,13 @@ 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);
 
+  // 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);
+
   save_grid_to_file(grid_file_name);
 }
 

nufi/save_results.h

@@ -2,12 +2,13 @@
 #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>
 
+class NuFISolver;
+
 void save_f(const NuFISolver &solver, unsigned int n,
             std::vector<GridStructure<1>> &grid_struct,
             std::vector<SolutionSnapshot<1>> &phi_history, unsigned int Nx_out,

src/nufi_solver.cc

@@ -4,6 +4,7 @@
 #include <boost/qvm/mat_access.hpp>
 #include <cmath>
 #include <cstdlib>
+#include <deal.II/base/exceptions.h>
 #include <deal.II/base/point.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/numerics/vector_tools.h>
@@ -49,6 +50,15 @@ std::vector<double> NuFISolver::eval_ftilda(
     for (size_t i = 0; i < x_size; ++i)
       X[i] = X[i] - Parameters::DT * U[i];
 
+    AssertThrow(
+        phi_history[n].solution.size() ==
+            grid_struct[phi_history[n].grid_version].dof_handler->n_dofs(),
+        ExcMessage("In eval_ftilda: Solution size = " +
+                   std::to_string(phi_history[n].solution.size()) +
+                   ", expected by grid_struct = " +
+                   std::to_string(grid_struct[phi_history[n].grid_version]
+                                      .dof_handler->n_dofs())));
+
     tmp = eval(X, grid_struct[phi_history[n].grid_version],
                phi_history[n].solution); // call eval only once