From e8d9fe9a40f0528e0254f59456f93f1b15a322dd Mon Sep 17 00:00:00 2001
From: "Vasco C. B. Ferreira" <vascocbf@gmail.com>
Date: Wed, 11 Mar 2026 14:05:39 +0100
Subject: [PATCH] added plotting fucntion for f

---
 nufi_solver.hpp     | 65 ++++++++++++++++++++++++++++++++++++++-------
 parameters.hpp      |  3 +++
 poisson_problem.hpp |  4 ---
 spline_field.hpp    |  2 --
 4 files changed, 59 insertions(+), 15 deletions(-)

diff --git a/nufi_solver.hpp b/nufi_solver.hpp
index f219da5..c8b7ae2 100644
--- a/nufi_solver.hpp
+++ b/nufi_solver.hpp
@@ -8,6 +8,7 @@
 #include <deal.II/numerics/fe_field_function.h>
 
 #include <iostream>
+#include <ostream>
 #include <vector>
 #include <cstddef>
 
@@ -25,13 +26,14 @@ public:
 
   void run();
   double eval_rho(unsigned int n, double x, const UniformSpline1D<double,4>& E_spline, unsigned int Nv = Parameters::NV);
+  double eval_ftilda(unsigned int n, double x, double u, const UniformSpline1D<double, 4>& E_spline);
+  void save_ftilda(unsigned int n, const UniformSpline1D<double,4>& E_spline, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename);
 
 private:
 
-  double eval_ftilda(unsigned int n, double x, double u, const UniformSpline1D<double, 4>& E_spline);
 
   unsigned int Nt = std::floor(Parameters::TMAX/Parameters::DT);
-  unsigned int Nx = Parameters::SPLINE_NX;
+  [[maybe_unused]] unsigned int Nx = Parameters::SPLINE_NX;
 
   double Lx = Parameters::LX;
 
@@ -112,10 +114,52 @@ private:
   const UniformSpline1D<double, 4> &E_spline;
 };
 
+inline void NuFISolver::save_ftilda(unsigned int n,
+                                    const UniformSpline1D<double,4>& E_spline,
+                                    unsigned int Nx_out,
+                                    unsigned int Nv_out,
+                                    const std::string &filename)
+{
+  std::ofstream file(filename);
+
+  double xmin = Parameters::X_DOMAIN_LEFT;
+  double xmax = Parameters::X_DOMAIN_RIGHT;
+
+  double vmin = Parameters::V_DOMAIN_LEFT;
+  double vmax = Parameters::V_DOMAIN_RIGHT;
+
+  double dx = (xmax - xmin) / Nx_out;
+  double dv = (vmax - vmin) / Nv_out;
+
+  file << Nx_out << " " << Nv_out << "\n";
+  file << xmin << " " << xmax << "\n";
+  file << vmin << " " << vmax << "\n";
+
+  for (unsigned int i = 0; i < Nx_out; ++i)
+  {
+      double x = xmin + (i + 0.5)*dx;
+
+      for (unsigned int j = 0; j < Nv_out; ++j)
+      {
+          double v = vmin + (j + 0.5)*dv;
+
+          double val = eval_ftilda(n, x, v, E_spline);
+
+          file << val;
+
+          if (j < Nv_out - 1)
+              file << " ";
+      }
+
+      file << "\n";
+  }
+
+  file.close();
+}
 
 inline void NuFISolver::run()
 {
-  std::cout << "Start of NuFISolver::run()\n";
+  std::cout << "Start of NuFISolver::run()\n\n";
   
   // init E_spline
 
@@ -129,17 +173,16 @@ inline void NuFISolver::run()
   for (unsigned int i=0; i<Nx; ++i)
   {
     [[maybe_unused]] double x = Parameters::X_DOMAIN_LEFT + i*dx;
-    E_grid[i] = 1; 
+    E_grid[i] = 0; 
   }
 
   UniformSpline1D<double,4> E_spline(E_grid, Parameters::X_DOMAIN_LEFT, Parameters::X_DOMAIN_RIGHT);
 
   for (unsigned int it = 0; it < Nt; ++it)
   {
-      std::cout << "Timestep " << it << " / " << Nt << std::endl;
+      std::cout << "Timestep " << it << " / " << Nt << std::endl << std::endl;
 
       // Step 1: Evaluate rho^n(x) using current E_spline
-      std::cout << "Start of eval_rho loop\n";
 
       rho.resize(Nx);
 
@@ -148,7 +191,6 @@ inline void NuFISolver::run()
           double x = (i + 0.5) * dx;
           rho[i] = eval_rho(it, x, E_spline, Parameters::NV);
       }
-      std::cout << "End of eval_rho loop\n";
 
       ChargeDensity_NuFI rho_function(*this, it, E_spline);
 
@@ -161,7 +203,12 @@ inline void NuFISolver::run()
 
       poisson.solve_step();
 
-      poisson.output_results(it);
+      if (it % Parameters::PLOT_FREQUENCY == 0)
+      {
+        std::cout << "Saving results... \n\n";
+        save_ftilda(it, E_spline, 128, 128, "results/ftilda_" + std::to_string(it) + ".dat");
+        poisson.output_results(it);
+      }
 
       E_grid = poisson.sample_electric_field(poisson, Nx, 0.0, Lx);
 
@@ -175,7 +222,7 @@ inline NuFISolver::NuFISolver()
   : order(Parameters::FE_DEGREE),
     poisson(order)
 {
-  std::cout << "Initializing Poisson\n";
+  std::cout << "Initializing dealii Poisson Solver\n";
   poisson.initialize();
 }
 #endif
diff --git a/parameters.hpp b/parameters.hpp
index ae43d0f..8fe7c69 100644
--- a/parameters.hpp
+++ b/parameters.hpp
@@ -29,6 +29,9 @@ namespace Parameters
 
   //spline options
   constexpr int SPLINE_NX = 256;
+
+  //Plotting options
+  constexpr int PLOT_FREQUENCY = 3;
 }
 
 #endif
diff --git a/poisson_problem.hpp b/poisson_problem.hpp
index 4566941..eb5a768 100644
--- a/poisson_problem.hpp
+++ b/poisson_problem.hpp
@@ -143,7 +143,6 @@ template<int dim>
 void PoissonProblem<dim>::create_mesh()
 {
 
-  std::cout << "Creating Mesh\n";
   GridGenerator::hyper_cube(triangulation,
                             Parameters::X_DOMAIN_LEFT,
                             Parameters::X_DOMAIN_RIGHT);
@@ -169,7 +168,6 @@ template <int dim>
 void PoissonProblem<dim>::setup_system()
 {
 
-  std::cout << "Setting up Poisson system\n";
   dof_handler.distribute_dofs(fe);
 
   constraints.clear();
@@ -222,7 +220,6 @@ public:
 template <int dim>
 void PoissonProblem<dim>::assemble_system()
 {
-  std::cout << "Assembling Poisson System\n";
   QGauss<dim>  quadrature_formula(fe.degree + 1);
   FEValues<dim> fe_values(fe, quadrature_formula,
                           update_values |
@@ -278,7 +275,6 @@ template <int dim>
 void PoissonProblem<dim>::solve()
 {
 
-  std::cout << "Calling PoissonProblem::solve()\n";
   SolverControl            solver_control(1000, 1e-12);
   SolverCG<Vector<double>> solver(solver_control);
 
diff --git a/spline_field.hpp b/spline_field.hpp
index fba9d05..cb9fff2 100644
--- a/spline_field.hpp
+++ b/spline_field.hpp
@@ -96,8 +96,6 @@ private:
 
     void interpolate(const std::vector<Real>& values)
     {
-
-      std::cout << "interpolating";
         size_t N = config.Nx;
 
         std::vector<Real> rhs(values);