nufi working, dealii solves poisson, E saved as independent spline

.gitignore

@@ -3,5 +3,6 @@ CMakeCache.txt
 compile_commands.json
 CMakeFiles/
 results/
+saved_sims/
 nufi_poisson
 *.ipynb

nufi/fields.h

@@ -47,19 +47,19 @@ double eval(double x, const double *coeffs) noexcept
     x -= Parameters::X_DOMAIN_LEFT;
 
     // Get "periodic position" in box at origin.
-    x = x - Parameters::LX * floor( x/Parameters::LX ); 
+    x = x - Parameters::LX * floor( x*Parameters::LX_INV ); 
 
     // Knot number
-    double x_knot = floor( x/Parameters::SPLINE_DX); 
+    double x_knot = floor( x*Parameters::SPLINE_DX_INV); 
 
     size_t ii = static_cast<size_t>(x_knot);
 
     // Convert x to reference coordinates.
-    x = x/Parameters::SPLINE_DX - x_knot;
+    x = x*Parameters::SPLINE_DX_INV - x_knot;
 
     // Scale according to derivative.
     double factor = 1;
-    for ( size_t i = 0; i < dx; ++i ) factor *= 1/Parameters::SPLINE_DX;
+    for ( size_t i = 0; i < dx; ++i ) factor *= 1*Parameters::SPLINE_DX_INV;
 
     return factor*splines1d::eval<double,Parameters::SPLINE_ORDER,dx>(x, coeffs + ii);
 }
@@ -142,7 +142,7 @@ void interpolate( real *coeffs, const real *values)
         coeffs[ i ] = tmp[ i % Parameters::SPLINE_NX ];
 }
 
-double integral_space_vector(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX)
+inline double integral_space_vector(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX)
 {
   double integral = 0.0;
   double x = Parameters::X_DOMAIN_LEFT;
@@ -153,7 +153,7 @@ double integral_space_vector(const double *current_coeffs, double dx = Parameter
   return integral;
 };
 
-double integral_space_vector_squared(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX)
+inline double integral_space_vector_squared(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX)
 {
   double integral = 0.0;
   double x = Parameters::X_DOMAIN_LEFT;

nufi/parameters.h

@@ -11,11 +11,13 @@ namespace Parameters
   constexpr double X_DOMAIN_LEFT  = 0.0;
   constexpr double X_DOMAIN_RIGHT = 4*M_PI;
   constexpr double LX = std::abs(X_DOMAIN_RIGHT- X_DOMAIN_LEFT);
+  constexpr double LX_INV = 1/LX;
 
   constexpr double V_DOMAIN_LEFT = -10.0;
   constexpr double V_DOMAIN_RIGHT = 10.0;
 
-  constexpr unsigned int NV = 562;
+  constexpr unsigned int NV = 512;
+  constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT)/NV;
 
   // deal.ii options
   constexpr unsigned int GLOBAL_REFINEMENT = 7;
@@ -25,15 +27,16 @@ namespace Parameters
 
   constexpr double EPS = 0.01;
   constexpr double WAVE_NR = 0.5;
-  constexpr double F0_FACTOR = 0.39894228040143267793994;
+  constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
 
   // NUFI options
-  constexpr double DT=1./16.;
+  constexpr double DT=1./10.;
-  constexpr unsigned int TMAX = 20;
+  constexpr unsigned int TMAX = 50;
 
   //spline options
-  constexpr int SPLINE_NX = 562;
+  constexpr int SPLINE_NX = 256;
-  constexpr double SPLINE_DX = LX/(SPLINE_NX-1);
+  constexpr double SPLINE_DX = LX/(SPLINE_NX);
+  constexpr double SPLINE_DX_INV = 1/SPLINE_DX;
   constexpr size_t SPLINE_ORDER = 4;
 
   //Plotting options

nufi/splines.h

@@ -1,5 +1,5 @@
-#ifndef SPLINES_H
+#ifndef SPLINES_HP
-#define SPLINES_H
+#define SPLINES_HP
 
 #include <cstddef>
 

src/nufi_solver.cc

@@ -27,8 +27,9 @@ double NuFISolver::eval_ftilda(unsigned int n,
 {
   if ( n == 0 ) return f0(x,u);
   
+  const size_t order = Parameters::SPLINE_ORDER;
   const size_t stride_x = 1;
-  const size_t stride_t = stride_x*(Nx + Parameters::SPLINE_ORDER - 1);
+  const size_t stride_t = stride_x*(Nx + order - 1);
 
   double Ex;
   const double *c;
@@ -107,26 +108,32 @@ void NuFISolver::run()
 
       std::vector<double> sampled_potential = poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of FE
 
-      save_space_vector(sampled_potential, "potential", it);
-
       // These have been tested to be equivalent
       // ////////////////////////////////////////////////
-      std::vector<double> E_vals = grad.compute(sampled_potential); // vector grad of FE solution
+      // std::vector<double> E_vals = grad.compute(sampled_potential); // vector grad of FE solution
-      save_space_vector(E_vals, "electric", it);
+      // save_space_vector(E_vals, "electric", it);
-      std::vector<double> E_vals_deal = poisson.sample_electric_field(x_min, x_max, Nx); // FE grad of soution
+      // std::vector<double> E_vals_deal = poisson.sample_electric_field(x_min, x_max, Nx); // FE grad of soution
-      save_space_vector(E_vals_deal, "electricdeal", it);
+      // save_space_vector(E_vals_deal, "electricdeal", it);
       // ////////////////////////////////////////////////                                             
 
-      double* current_coeffs = coeffs.get() + it*stride_t;
-      interpolate<double, Parameters::SPLINE_ORDER>(current_coeffs, E_vals.data());
 
+      // interpolate and save current field
+      double* current_coeffs = coeffs.get() + it*stride_t;
+      interpolate<double, Parameters::SPLINE_ORDER>(current_coeffs, sampled_potential.data());
+
+      std::vector<double> E_x(Nx,0.0) ;
+      for(size_t ix=0; ix<Nx; ++ix)
+      {
+        E_x[ix] = -eval<1>(Parameters::X_DOMAIN_LEFT+ix*dx, current_coeffs);
+      }
 
       if (it % Parameters::PLOT_FREQUENCY == 0)
       {
         std::cout << "Saving results... \n\n";
         save_ftilda(*this, it, coeffs.get(), 128, 128, "results/ftilda_" + std::to_string(it) + ".dat");
         save_rho(*this, it, coeffs.get(), 128, "results/rho_" + std::to_string(it) + ".dat");
-        save_Efield(it, coeffs.get(), 128, "results/field_" + std::to_string(it) + ".dat");
+        // save_Efield(it, coeffs.get(), 128, "results/field_" + std::to_string(it) + ".dat");
+        save_space_vector(E_x, "field", it);
 
         double int_val = 0.5 * integral_space_vector_squared(current_coeffs);
         int_E_squared.push_back(int_val);