initial working 1x2v working, review and optimization needed

Makefile

@@ -67,10 +67,10 @@ RM = /usr/local/bin/cmake -E rm -f
 EQUALS = =
 
 # The top-level source directory on which CMake was run.
-CMAKE_SOURCE_DIR = /home/vf/Thesis/dealII_NuFI/1x1v
+CMAKE_SOURCE_DIR = /home/vf/Thesis/dealII_NuFI
 
 # The top-level build directory on which CMake was run.
-CMAKE_BINARY_DIR = /home/vf/Thesis/dealII_NuFI/1x1v
+CMAKE_BINARY_DIR = /home/vf/Thesis/dealII_NuFI
 
 #=============================================================================
 # Targets provided globally by CMake.
@@ -89,7 +89,7 @@ rebuild_cache/fast: rebuild_cache
 # Special rule for the target edit_cache
 edit_cache:
 	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake cache editor..."
-	/usr/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
+	/usr/local/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
 .PHONY : edit_cache
 
 # Special rule for the target edit_cache
@@ -99,9 +99,9 @@ edit_cache/fast: edit_cache
 
 # The main all target
 all: cmake_check_build_system
-	$(CMAKE_COMMAND) -E cmake_progress_start /home/vf/Thesis/dealII_NuFI/1x1v/CMakeFiles /home/vf/Thesis/dealII_NuFI/1x1v//CMakeFiles/progress.marks
+	$(CMAKE_COMMAND) -E cmake_progress_start /home/vf/Thesis/dealII_NuFI/CMakeFiles /home/vf/Thesis/dealII_NuFI//CMakeFiles/progress.marks
 	$(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 all
-	$(CMAKE_COMMAND) -E cmake_progress_start /home/vf/Thesis/dealII_NuFI/1x1v/CMakeFiles 0
+	$(CMAKE_COMMAND) -E cmake_progress_start /home/vf/Thesis/dealII_NuFI/CMakeFiles 0
 .PHONY : all
 
 # The main clean target

README.md

@@ -5,9 +5,9 @@ This simulation of the Vlasov-Poisson system dimensions uses
 - [deal.ii](https://dealii.org/) FEM package
 
 ---
-dimensions: 1x1v
+dimensions: 1x2v
 
-status:
+status: working, needs testing, not optimized (at all!)
 
 Refinement working:
 - grid versions saved on a vector

cmake_install.cmake

@@ -1,4 +1,4 @@
-# Install script for directory: /home/vf/Thesis/dealII_NuFI/1x1v
+# Install script for directory: /home/vf/Thesis/dealII_NuFI
 
 # Set the install prefix
 if(NOT DEFINED CMAKE_INSTALL_PREFIX)
@@ -50,5 +50,5 @@ endif()
 
 string(REPLACE ";" "\n" CMAKE_INSTALL_MANIFEST_CONTENT
        "${CMAKE_INSTALL_MANIFEST_FILES}")
-file(WRITE "/home/vf/Thesis/dealII_NuFI/1x1v/${CMAKE_INSTALL_MANIFEST}"
+file(WRITE "/home/vf/Thesis/dealII_NuFI/${CMAKE_INSTALL_MANIFEST}"
      "${CMAKE_INSTALL_MANIFEST_CONTENT}")

nufi/fields.h

@@ -1,7 +1,7 @@
 #ifndef FIELDS_H
 #define FIELDS_H
 
-#include "grids.h"
+#include "nufi/grids.h"
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"
 #include <cmath>
@@ -21,26 +21,15 @@ inline std::vector<double> make_x_eval(size_t Nx) {
   for (unsigned int i = 0; i < Nx; ++i) {
     x_eval_E.push_back(Parameters::X_DOMAIN_LEFT + (i + 0.5) * dx);
   }
-  // std::vector<double> x_eval(Nx);
-  // const double dx = Parameters::LX / Nx;
-  // for (size_t i = 0; i < Nx; ++i)
-  //   x_eval[i] = Parameters::X_DOMAIN_LEFT + i * dx;
   return x_eval_E;
 }
 
-inline void reset_x_eval(std::vector<double> &x_vals) {
+inline double f0(const double x, const double v1, const double v2,
-  const size_t Nx = x_vals.size();
-  const double dx = Parameters::LX / Nx;
-  for (size_t i = 0; i < Nx; ++i)
-    x_vals[i] = Parameters::X_DOMAIN_LEFT + i * dx;
-};
-
-inline double f0(const double x, const double v,
                  const double eps = Parameters::EPS,
                  const double k = Parameters::WAVE_NR) {
   const double prefactor =
       Parameters::F0_FACTOR * (1.0 + eps * std::cos(k * x));
-  const double gaussian = v * v * std::exp(-0.5 * v * v);
+  const double gaussian = v1 * v1 * std::exp(-0.5 * (v1 * v1 + v2 * v2));
 
   return prefactor * gaussian;
 }

nufi/nufi_solver.h

@@ -24,13 +24,14 @@ public:
   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;
+           const unsigned int Nv1 = Parameters::NV_1,
+           const unsigned int Nv2 = Parameters::NV_2) const;
   std::vector<double>
-  eval_ftilda(unsigned int n, std::vector<double> x, double u,
+  eval_ftilda(unsigned int n, std::vector<double> x, double u1, double u2,
               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,
+  eval_f(unsigned int n, std::vector<double> x, double u1, double u2,
          const std::vector<GridStructure<1>> &grid_struct,
          const std::vector<SolutionSnapshot<1>> &phi_history) const;
 
@@ -38,7 +39,7 @@ public:
   eval_rho_points(unsigned int n, const std::vector<Point<1>> &points,
                   const std::vector<GridStructure<1>> &grid_struct,
                   const std::vector<SolutionSnapshot<1>> &phi_history,
-                  const unsigned int Nv) const;
+                  const unsigned int Nv1, const unsigned int Nv2) const;
 
 private:
   unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT);

nufi/parameters.h

@@ -16,12 +16,17 @@ constexpr double LX_INV = 1 / LX;
 constexpr size_t CALC_NX = 256;
 constexpr double CALC_DX = LX / CALC_NX;
 
-constexpr double V_DOMAIN_LEFT = -10.;
+constexpr double V_DOMAIN_LEFT_1 = -10.;
-constexpr double V_DOMAIN_RIGHT = 10.;
+constexpr double V_DOMAIN_RIGHT_1 = 10.;
 
-constexpr unsigned int NV = 128;
+constexpr unsigned int NV_1 = 128;
-constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV;
+constexpr double DV_1 = std::abs(V_DOMAIN_RIGHT_1 - V_DOMAIN_LEFT_1) / NV_1;
 
+constexpr double V_DOMAIN_LEFT_2 = -10.;
+constexpr double V_DOMAIN_RIGHT_2 = 10.;
+
+constexpr unsigned int NV_2 = 128;
+constexpr double DV_2 = std::abs(V_DOMAIN_RIGHT_2 - V_DOMAIN_LEFT_2) / NV_2;
 // deal.ii options
 constexpr unsigned int GLOBAL_REFINEMENT = 6;
 constexpr unsigned int FE_DEGREE = 3;
@@ -45,6 +50,7 @@ constexpr unsigned int REFINE_FREQUENCY = 30;
 constexpr int PLOT_FREQUENCY = 10;
 constexpr size_t PLOT_NX = CALC_NX;
 constexpr double PLOT_DX = LX / PLOT_NX;
+constexpr double PLOT_FIXED_V2 = 0.;
 const std::string PLOT_DIR = "results/";
 } // namespace Parameters
 

nufi/save_results.h

@@ -2,14 +2,13 @@
 #define SAVE_RESULTS_H
 
 #include "nufi/grids.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,
+void save_f(const NuFISolver &solver, unsigned int n, const double v2_0,
             std::vector<GridStructure<1>> &grid_struct,
             std::vector<SolutionSnapshot<1>> &phi_history, unsigned int Nx_out,
             unsigned int Nv_out, const std::string &filename);

src/nufi_solver.cc

@@ -12,7 +12,6 @@
 #include <cstddef>
 #include <fstream>
 #include <iostream>
-#include <memory>
 #include <ostream>
 #include <string>
 #include <vector>
@@ -27,18 +26,18 @@
 using namespace dealii;
 
 std::vector<double> NuFISolver::eval_ftilda(
-    unsigned int n, std::vector<double> X, double u,
+    unsigned int n, std::vector<double> X, double u1, double u2,
     const std::vector<GridStructure<1>> &grid_struct,
     const std::vector<SolutionSnapshot<1>> &phi_history) const {
 
   size_t x_size = X.size();
 
-  std::vector<double> U(x_size, u);
+  std::vector<double> U1(x_size, u1);
+  std::vector<double> U2(x_size, u2);
   std::vector<double> results(x_size);
   if (n == 0) {
     for (size_t i = 0; i < x_size; ++i)
-      results[i] = f0(X[i], U[i]);
+      results[i] = f0(X[i], U1[i], U2[i]);
-    // reset_x_eval(X);
     return results;
   }
 
@@ -48,61 +47,53 @@ std::vector<double> NuFISolver::eval_ftilda(
   // We omit the initial half-step.
   while (--n) {
     for (size_t i = 0; i < x_size; ++i)
-      X[i] = X[i] - Parameters::DT * U[i];
+      X[i] = X[i] - Parameters::DT * U1[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())));
-    AssertThrow(
-        grid_struct[phi_history[n].grid_version].grid_version ==
-            phi_history[n].grid_version,
-        ExcMessage(
-            "grid.grid_version not equal to phi_history[n].grid_version"));
 
     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];
+      U1[i] = U1[i] + Parameters::DT * Ex[i];
+
+      // IMPORTANT!! LOOK at this before 1x2v, and check with Paul
+      //
+      // No change to U2 because dv1/dt = E_x
+      // but dv2/dt = E_y = 0 (because were in 1x not 2x)
+      // U2[i] = U2[i] + Parameters::DT * Ey[i]; TO_BUILD
     }
   }
 
   // The final half-step.
   for (size_t i = 0; i < x_size; ++i)
-    X[i] = X[i] - Parameters::DT * U[i];
+    X[i] = X[i] - Parameters::DT * U1[i];
 
   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];
+    U1[i] = U1[i] + 0.5 * Parameters::DT * Ex[i];
   }
   for (size_t i = 0; i < x_size; ++i)
-    results[i] = f0(X[i], U[i]);
+    results[i] = f0(X[i], U1[i], U2[i]);
-  // reset_x_eval(X);
   return results;
 }
 
 std::vector<double>
-NuFISolver::eval_f(unsigned int n, std::vector<double> X, double u,
+NuFISolver::eval_f(unsigned int n, std::vector<double> X, double u1, double u2,
                    const std::vector<GridStructure<1>> &grid_struct,
                    const std::vector<SolutionSnapshot<1>> &phi_history) const {
 
   size_t x_size = X.size();
 
-  std::vector<double> U(x_size, u);
+  std::vector<double> U1(x_size, u1);
+  std::vector<double> U2(x_size, u2);
+
   std::vector<double> results(x_size);
   if (n == 0) {
     for (size_t i = 0; i < x_size; ++i)
-      results[i] = f0(X[i], U[i]);
+      results[i] = f0(X[i], U1[i], U2[i]);
-    // reset_x_eval(X);
     return results;
   }
 
@@ -114,35 +105,34 @@ NuFISolver::eval_f(unsigned int n, std::vector<double> X, double u,
              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];
+    U1[i] = U1[i] + 0.5 * Parameters::DT * Ex[i];
   }
 
   while (--n) {
     for (size_t i = 0; i < x_size; ++i)
-      X[i] = X[i] - Parameters::DT * U[i];
+      X[i] = X[i] - Parameters::DT * U1[i];
 
     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];
+      U1[i] = U1[i] + Parameters::DT * Ex[i];
     }
   }
 
   // The final half-step.
   for (size_t i = 0; i < x_size; ++i)
-    X[i] = X[i] - Parameters::DT * U[i];
+    X[i] = X[i] - Parameters::DT * U1[i];
 
   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];
+    U1[i] = U1[i] + 0.5 * Parameters::DT * Ex[i];
   }
 
   for (size_t i = 0; i < x_size; ++i)
-    results[i] = f0(X[i], U[i]);
+    results[i] = f0(X[i], U1[i], U2[i]);
-  // reset_x_eval(X);
   return results;
 }
 
@@ -150,26 +140,30 @@ std::vector<double>
 NuFISolver::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) const {
+                     const unsigned int Nv1, const unsigned int Nv2) const {
   size_t x_size = X.size();
 
-  const double dv =
+  const double dv_0 = Parameters::DV_1;
-      (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
+  const double v_min_0 = Parameters::V_DOMAIN_LEFT_1 + 0.5 * dv_0;
 
-  const double v_min = Parameters::V_DOMAIN_LEFT + 0.5 * dv;
+  const double dv_1 = Parameters::DV_2;
+  const double v_min_1 = Parameters::V_DOMAIN_LEFT_2 + 0.5 * dv_1;
 
   std::vector<double> integral(x_size, 0.0);
   std::vector<double> tmp_int(x_size);
 
-  for (unsigned int i = 0; i < Nv; ++i) {
+  for (unsigned int i = 0; i < Nv1; ++i) {
-    tmp_int = eval_ftilda(n, X, v_min + i * dv, grid_struct,
+    for (unsigned int j = 0; j < Nv2; ++j) {
-                          phi_history); // used eval_ftilda once per i
+      tmp_int =
-    for (size_t ii = 0; ii < x_size; ++ii)
+          eval_ftilda(n, X, v_min_0 + i * dv_0, v_min_1 + j * dv_1, grid_struct,
-      integral[ii] += tmp_int[ii];
+                      phi_history); // used eval_ftilda once per i
+      for (size_t ii = 0; ii < x_size; ++ii)
+        integral[ii] += tmp_int[ii];
+    }
   }
 
   for (size_t i = 0; i < x_size; ++i)
-    integral[i] = 1 - integral[i] * dv;
+    integral[i] = 1 - integral[i] * dv_0;
   return integral;
 }
 
@@ -177,9 +171,11 @@ std::vector<double>
 NuFISolver::eval_rho_points(unsigned int n, const std::vector<Point<1>> &points,
                             const std::vector<GridStructure<1>> &grid_struct,
                             const std::vector<SolutionSnapshot<1>> &phi_history,
-                            const unsigned int Nv) const {
+                            const unsigned int Nv1,
+                            const unsigned int Nv2) const {
   std::vector<double> point_vector = Point_vector_to_double_vector(points);
-  return NuFISolver::eval_rho(n, point_vector, grid_struct, phi_history, Nv);
+  return NuFISolver::eval_rho(n, point_vector, grid_struct, phi_history, Nv1,
+                              Nv2);
 }
 
 void NuFISolver::run() {
@@ -280,7 +276,8 @@ void NuFISolver::run() {
       for (size_t i = 0; i < points.size(); ++i)
         x[i] = points[i][0];
 
-      return eval_rho(it, x, grid_versions, phi_history, Parameters::NV);
+      return eval_rho(it, x, grid_versions, phi_history, Parameters::NV_1,
+                      Parameters::NV_2);
     });
 
     if (it % Parameters::REFINE_FREQUENCY == 0) {
@@ -306,8 +303,9 @@ void NuFISolver::run() {
       double plot_start = timer.elapsed();
 
       std::cout << "Saving results...   ";
-      save_f(*this, it, grid_versions, phi_history, Parameters::PLOT_NX,
+      save_f(*this, it, Parameters::PLOT_FIXED_V2, grid_versions, phi_history,
-             Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
+             Parameters::PLOT_NX, Parameters::NV_1,
+             "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");
 

src/save_results.cc

@@ -12,7 +12,7 @@
 #include <string>
 #include <vector>
 
-void save_f(const NuFISolver &solver, unsigned int n,
+void save_f(const NuFISolver &solver, unsigned int n, const double v2_0,
             std::vector<GridStructure<1>> &grid_struct,
             std::vector<SolutionSnapshot<1>> &phi_history, unsigned int Nx_out,
             unsigned int Nv_out, const std::string &filename) {
@@ -21,21 +21,21 @@ void save_f(const NuFISolver &solver, unsigned int n,
   double xmin = Parameters::X_DOMAIN_LEFT;
   double xmax = Parameters::X_DOMAIN_RIGHT;
 
-  double vmin = Parameters::V_DOMAIN_LEFT;
+  double vmin_1 = Parameters::V_DOMAIN_LEFT_1;
-  double vmax = Parameters::V_DOMAIN_RIGHT;
+  double vmax_1 = Parameters::V_DOMAIN_RIGHT_1;
 
-  double dv = (vmax - vmin) / Nv_out;
+  double dv_1 = (vmax_1 - vmin_1) / Nv_out;
 
   file << Nx_out << " " << Nv_out << "\n";
   file << xmin << " " << xmax << "\n";
-  file << vmin << " " << vmax << "\n";
+  file << vmin_1 << " " << vmax_1 << "\n";
 
   std::vector<double> x_eval = make_x_eval(Nx_out);
   std::vector<double> val(Nx_out);
 
   for (unsigned int j = 0; j < Nv_out; ++j) {
-    double v = vmin + (j + 0.5) * dv;
+    double v = vmin_1 + (j + 0.5) * dv_1;
-    val = solver.eval_f(n, x_eval, v, grid_struct, phi_history);
+    val = solver.eval_f(n, x_eval, v, v2_0, grid_struct, phi_history);
 
     for (unsigned int i = 0; i < Nx_out; ++i) {
       file << val[i];