added std::vector gradient

README.md

@@ -6,4 +6,4 @@ This simulation of the Vlasov-Poisson system in 1x1v dimensions uses
 
 --- 
 Todo:
-- Correct E field sampler from dealii's sovler
+- Correct spline interpolator/evaluation for discontinuity

nufi/fields.h

@@ -142,6 +142,41 @@ void interpolate( real *coeffs, const real *values)
         coeffs[ i ] = tmp[ i % Parameters::SPLINE_NX ];
 }
 
+class Gradient {
+public:
+    Gradient(double xmin, double xmax, unsigned int Nx)
+        : xmin_(xmin), xmax_(xmax), Nx_(Nx)
+    {
+        if (xmax_ <= xmin_) {
+            throw std::invalid_argument("xmax must be greater than xmin");
+        }
+    }
+
+    std::vector<double> compute(const std::vector<double>& values) const {
+        size_t n = values.size();
+        if (n < 2) {
+            throw std::invalid_argument("Need at least 2 points");
+        }
+
+        std::vector<double> grad(n);
+
+        double dx = (xmax_ - xmin_) / (n-1);
+        // periodic boundaries
+        grad[0]     = -(values[1]   - values[n-1]) / (2.0 * dx);
+        grad[n-1]   = -(values[0]   - values[n-2]) / (2.0 * dx);
+
+        for (size_t i = 1; i < n-1; ++i) {
+            grad[i] = -(values[i+1] - values[i-1]) / (2.0 * dx);
+        }
+
+
+        return grad;
+    }
+private:
+    double xmin_;
+    double xmax_;
+    [[maybe_unused]] unsigned int Nx_;
+};
 
 template <int dim>
 class ChargeDensity : public Function<dim> // only uses f0

nufi/nufi_solver.h

@@ -1,6 +1,7 @@
 #ifndef NUFI_SOLVER_H
 #define NUFI_SOLVER_H
 
+#include <boost/qvm/mat_access.hpp>
 #include <vector>
 #include <cmath>
 #include <deal.II/base/point.h>
@@ -30,6 +31,9 @@ private:
 
   double Lx = Parameters::LX;
 
+  double x_min = Parameters::X_DOMAIN_LEFT;
+  double x_max = Parameters::X_DOMAIN_RIGHT;
+
   std::vector<double> rho;
 
   unsigned int order;

nufi/poisson_problem.h

@@ -3,6 +3,8 @@
 
 #include <deal.II/base/function.h>
 
+#include <deal.II/base/mpi_remote_point_evaluation.h>
+#include <deal.II/base/point.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/template_constraints.h>
@@ -34,6 +36,7 @@
 #include <deal.II/numerics/matrix_tools.h>
 #include <deal.II/numerics/fe_field_function.h>
 
+#include <deal.II/numerics/vector_tools_evaluate.h>
 #include <deal.II/numerics/vector_tools_interpolate.h>
 #include <memory>
 #include <string>
@@ -61,9 +64,8 @@ public:
   const Vector<double> &get_solution() const { return solution; }
   const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
 
-  std::vector<double> sample_electric_field(unsigned int Nx,
+  std::vector<double> sample_electric_field(double x_min, double x_max, unsigned int Nx);
-                                            double x_min,
+  std::vector<double> sample_electric_potential(double x_min, double x_max, unsigned int Nx);
-                                            double x_max);
 
 private:
   void create_mesh();
@@ -104,33 +106,71 @@ PoissonProblem<dim>::PoissonProblem(unsigned int degree)
 {}
 
 template <int dim>
-std::vector<double> PoissonProblem<dim>::sample_electric_field(
+std::vector<double> PoissonProblem<dim>::sample_electric_field(double x_min,double x_max,unsigned int Nx)
-    unsigned int Nx,
-    double x_min,
-    double x_max)
 {
-  this -> get_solution();
+  std::vector<double> E_values(Nx);
-  this -> get_dof_handler();
 
-  Functions::FEFieldFunction<dim, Vector<double>>
+  const double dx = (x_max - x_min) / (Nx - 1);
-    field_function(dof_handler, solution, mapping);
 
+  for (unsigned int i = 0; i < Nx; ++i)
+  {
+      const double x = x_min + i * dx;
+      const Point<dim> point(x);
+
+      // 1. Find the active cell containing x
+      const auto cell_point_pair =
+          GridTools::find_active_cell_around_point(mapping,
+                                                   dof_handler,
+                                                   point);
+
+      const auto cell = cell_point_pair.first;
+      const Point<dim> &unit_point = cell_point_pair.second;
+
+      // 2. FEPointEvaluation expects an ArrayView of points
+      std::vector<Point<dim>> points(1, unit_point);
+      ArrayView<const Point<dim>> point_view(points);
+
+      FEPointEvaluation<1, dim> evaluator(mapping,
+                                          dof_handler.get_fe(),
+                                          update_gradients);
+
+      // reinit with ArrayView of points
+      evaluator.reinit(cell, point_view);
+
+      Vector<double> local_dofs(dof_handler.get_fe().dofs_per_cell);
+      cell->get_dof_values(solution, local_dofs);
+
+      // 3. Evaluate gradient at this point
+      evaluator.evaluate(local_dofs, EvaluationFlags::gradients);
+
+      const Tensor<1, dim> grad_phi = evaluator.get_gradient(0);
+
+      // 4. Compute E = -grad(phi)
+      E_values[i] = -grad_phi[0];
+  }
+
+  return E_values;
+}
+
+template <int dim>
+std::vector<double> PoissonProblem<dim>::sample_electric_potential(
+    double x_min,
+    double x_max,
+    unsigned int Nx)
+{
   std::vector<double> values(Nx);
+  std::vector<Point<dim>> eval_points(Nx);
 
   double Lx = x_max - x_min;
   double dx = Lx / Nx;
 
-  for (unsigned int i = 0; i < Nx; ++i)
+  for(unsigned int i=0 ; i<Nx; ++i)
-  {
+    eval_points[i] = Point<1, double>(x_min + i * dx);
-      double x = x_min + i * dx;
 
-      Point<dim> p;
+  Utilities::MPI::RemotePointEvaluation<dim,dim> cache;
-      p[0] = x;
+  cache.reinit(eval_points, triangulation, mapping);
 
-      Tensor<1, dim> grad = field_function.gradient(p);
+  values = VectorTools::point_values<dim>(cache, dof_handler, solution);
-
-      values[i] = -grad[0];  // E = -dφ/dx
-  }
 
   return values;
 }

src/nufi_solver.cc

@@ -28,7 +28,7 @@ double NuFISolver::eval_ftilda(unsigned int n,
   if ( n == 0 ) return f0(x,u);
 
   const size_t stride_x = 1;
-  const size_t stride_t = stride_x*(Parameters::SPLINE_NX + Parameters::SPLINE_ORDER - 1);
+  const size_t stride_t = stride_x*(Nx + Parameters::SPLINE_ORDER - 1);
 
   double Ex;
   const double *c;
@@ -81,6 +81,8 @@ void NuFISolver::run()
 
   if ( rho == nullptr ) throw std::bad_alloc {};
 
+  Gradient grad(x_min, x_max, Nx);
+
   for (unsigned int it = 0; it < Nt; ++it)
   {
       std::cout << "Timestep " << it << " / " << Nt << std::endl << std::endl;
@@ -97,10 +99,17 @@ void NuFISolver::run()
     		rho.get()[i] =  ith_rho;
     	}
 
-      poisson.set_rhs_function(std::make_unique<ChargeDensity_NuFI<1>>(rho.get(), Parameters::SPLINE_NX));
+      poisson.set_rhs_function(std::make_unique<ChargeDensity_NuFI<1>>(rho.get(), Nx));
       poisson.solve_step();
 
-      std::vector<double> E_vals = poisson.sample_electric_field(Parameters::SPLINE_NX, Parameters::X_DOMAIN_LEFT, Parameters::X_DOMAIN_RIGHT);
+      std::vector<double> sampled_potential = poisson.sample_electric_potential(x_min, x_max, Nx);
+
+      save_space_vector(sampled_potential, "potential", it);
+
+      std::vector<double> E_vals = grad.compute(sampled_potential); 
+      // std::vector<double> E_vals = poisson.sample_electric_field(x_min, x_max, Nx);
+      
+      save_space_vector(E_vals, "electric", it);
 
       double* current_coeffs = coeffs.get() + it*stride_t;
       interpolate<double, Parameters::SPLINE_ORDER>(current_coeffs, E_vals.data());

src/save_results.cc

@@ -2,7 +2,9 @@
 
 #include "nufi/parameters.h"
 #include <fstream>
+#include <stdexcept>
 #include <string>
+#include <vector>
 #include "nufi/nufi_solver.h"