From bf7ad1e59513f974db757c4fc091725962607e41 Mon Sep 17 00:00:00 2001
From: "Vasco C. B. Ferreira" <vascocbf@gmail.com>
Date: Sun, 8 Mar 2026 17:05:18 +0100
Subject: [PATCH] nufi working but slow, need to save field from previous time
 step and/or parallel loop to solve eval_rho

---
 nufi_solver.hpp     | 11 +++++++----
 poisson_problem.hpp |  9 ++++++++-
 2 files changed, 15 insertions(+), 5 deletions(-)

diff --git a/nufi_solver.hpp b/nufi_solver.hpp
index 3791510..f8b7625 100644
--- a/nufi_solver.hpp
+++ b/nufi_solver.hpp
@@ -91,10 +91,12 @@ inline double NuFISolver::eval_ftilda(unsigned int n,
  
   double Ex;
 
+
   // Initial half-step.
   Ex = evaluate_E(x);
   u += 0.5*dt*Ex;
 
+
   while ( --n )
   {
       x -= dt*u;
@@ -106,7 +108,6 @@ inline double NuFISolver::eval_ftilda(unsigned int n,
   x -= dt*u;
   Ex = evaluate_E(x);
   u += 0.5*dt*Ex; // is this line useless ?
-  
   double x_periodic = x - Lx * std::floor(x / Lx);
   double u_periodic = u - Lu * std::floor(u / Lu);
 
@@ -120,7 +121,6 @@ inline double NuFISolver::eval_rho(const unsigned int n,
   const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
 
   double integral = 0.0;
-
   for (unsigned int i = 0; i < Nv; ++i)
   {
     const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv;
@@ -170,13 +170,15 @@ inline void NuFISolver::run()
 
 
     double dx = Lx / Nx;
-
+    
+    std::cout << "Start of eval_rho step with Nx = "<< Nx<< "\n";
     for (unsigned int i = 0; i < Nx; ++i)
     {
       double x = (i + 0.5) * dx;
       rho[i] = eval_rho(n, x);
     }
-
+    std::cout << "End of eval_rho step\n";
+        
     solve_poisson(n);
   }
 
@@ -187,6 +189,7 @@ inline NuFISolver::NuFISolver()
   : order(Parameters::FE_DEGREE),
     poisson(order, Parameters::NV)
 {
+  std::cout << "Initializing Poisson\n";
   poisson.initialize();
 
   Nx = poisson.get_dof_handler().n_dofs();
diff --git a/poisson_problem.hpp b/poisson_problem.hpp
index b5b6f57..472a987 100644
--- a/poisson_problem.hpp
+++ b/poisson_problem.hpp
@@ -108,6 +108,8 @@ PoissonProblem<dim>::PoissonProblem(unsigned int degree, unsigned int Nv)
 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);
@@ -132,6 +134,8 @@ void PoissonProblem<dim>::create_mesh()
 template <int dim>
 void PoissonProblem<dim>::setup_system()
 {
+
+  std::cout << "Setting up Poisson system\n";
   dof_handler.distribute_dofs(fe);
 
   constraints.clear();
@@ -184,6 +188,7 @@ 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 |
@@ -238,6 +243,8 @@ void PoissonProblem<dim>::assemble_system()
 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);
 
@@ -313,7 +320,7 @@ void PoissonProblem<dim>::solve_step()
 {
   system_matrix = 0;
   system_rhs = 0;
-
+  std::cout << "Calling PoissonProblem::solve_step()\n";
   assemble_system();
   solve();
 }