refinment working with history of grids

nufi/grids.h

@@ -48,6 +48,8 @@ template <int dim> struct GridStructure {
 
     std::vector<double> values(points.size());
 
+    // Uses point_gradient
+
     // #pragma omp parallel for
     //     for (unsigned int p = 0; p < points.size(); ++p) {
     //
@@ -60,6 +62,7 @@ template <int dim> struct GridStructure {
     //     return values;
     //   }
 
+    // Uses cell locator
 #pragma omp parallel
     {
       std::vector<double> local_solution_buffer(fe->n_dofs_per_cell());
@@ -82,9 +85,6 @@ template <int dim> struct GridStructure {
         values[p] = evaluator.get_gradient(0)[0];
       }
     }
-    // AssertThrow(dof_handler->n_dofs() == solution.size(),
-    //             ExcMessage("Solution vector size does not match
-    //             DoFHandler.)");
     return values;
   }
 };
@@ -102,21 +102,11 @@ GridStructure<dim> make_grid_snapshot(PoissonProblem<dim> &poisson) {
   grid.fe = std::make_unique<FE_Q<dim>>(poisson.get_fe());
 
   grid.dof_handler = std::make_unique<DoFHandler<dim>>(*grid.triangulation);
-
   grid.dof_handler->distribute_dofs(*grid.fe);
 
   grid.mapping = std::make_unique<MappingQ<dim>>(poisson.get_mapping());
-  grid.constraints =
-      std::make_unique<AffineConstraints<double>>(poisson.get_constraints());
 
-  grid.locator = std::make_unique<CellLocator<dim>>();
-
-  // START: transfer poisson.solution to saved grid dofs
-  SolutionTransfer<dim> solution_transfer(*grid.dof_handler);
-  const Vector<double> coarse_solution = poisson.solution;
-  solution_transfer.prepare_for_coarsening_and_refinement(coarse_solution);
-  // START: setup_system();
-  grid.constraints->clear();
+  grid.constraints = std::make_unique<AffineConstraints<double>>();
 
   DoFTools::make_hanging_node_constraints(*grid.dof_handler, *grid.constraints);
   DoFTools::make_periodicity_constraints(*grid.dof_handler, 0, 1, 0,
@@ -147,21 +137,10 @@ GridStructure<dim> make_grid_snapshot(PoissonProblem<dim> &poisson) {
 
   grid.constraints->add_line(gauge_dof);
   grid.constraints->set_inhomogeneity(gauge_dof, 0.0);
-
   grid.constraints->close();
 
-  DynamicSparsityPattern dsp(grid.dof_handler->n_dofs());
-  DoFTools::make_sparsity_pattern(*grid.dof_handler, dsp, *grid.constraints);
-  poisson.sparsity_pattern.copy_from(dsp);
-  poisson.system_matrix.reinit(poisson.sparsity_pattern);
-  poisson.solution.reinit(grid.dof_handler->n_dofs());
-  poisson.system_rhs.reinit(grid.dof_handler->n_dofs());
-
+  grid.locator = std::make_unique<CellLocator<dim>>();
   grid.locator->rebuild(*grid.dof_handler, *grid.triangulation);
-  // END
-
-  solution_transfer.interpolate(coarse_solution, poisson.solution);
-  // END
 
   // START: diagnostics
   AssertThrow(
@@ -169,32 +148,18 @@ GridStructure<dim> make_grid_snapshot(PoissonProblem<dim> &poisson) {
           grid.constraints->n_constraints(),
       ExcMessage(
           "PoissonProblem constraints doesn't match Snapshot constraints"));
-  for (auto c1 = poisson.get_dof_handler().begin_active(),
-            c2 = grid.dof_handler->begin_active();
-       c1 != poisson.get_dof_handler().end(); ++c1, ++c2) {
-    std::vector<types::global_dof_index> d1(c1->get_fe().dofs_per_cell);
-    std::vector<types::global_dof_index> d2(c2->get_fe().dofs_per_cell);
-
-    c1->get_dof_indices(d1);
-    c2->get_dof_indices(d2);
-
-    AssertThrow(d1 == d2, ExcInternalError());
-  }
+  // for (auto c1 = poisson.get_dof_handler().begin_active(),
+  //           c2 = grid.dof_handler->begin_active();
+  //      c1 != poisson.get_dof_handler().end(); ++c1, ++c2) {
+  //   std::vector<types::global_dof_index> d1(c1->get_fe().dofs_per_cell);
+  //   std::vector<types::global_dof_index> d2(c2->get_fe().dofs_per_cell);
   //
-  // auto support_points =
-  //     DoFTools::map_dofs_to_support_points(*grid.mapping, *grid.dof_handler);
+  //   c1->get_dof_indices(d1);
+  //   c2->get_dof_indices(d2);
   //
-  // std::cout << "COPY\n";
-  //
-  // for (const auto &[dof, point] : support_points) {
-  //   std::cout << dof << " : " << point[0] << "\n";
+  //   AssertThrow(d1 == d2, ExcInternalError());
   // }
-  // for (auto cell : grid.dof_handler->active_cell_iterators()) {
-  //   std::cout << "@ GridStructure: " << cell->id() << "  " <<
-  //   cell->center()[0]
-  //             << '\n';
-  // }
-  // END
+  // END: diagnostics
 
   return grid;
 }

nufi/nufi_solver.h

@@ -26,11 +26,11 @@ public:
            const std::vector<SolutionSnapshot<1>> &phi_history,
            const unsigned int Nv = Parameters::NV) const;
   std::vector<double>
-  eval_ftilda(unsigned int, std::vector<double> &x, double u,
+  eval_ftilda(unsigned int n, std::vector<double> x, double u,
               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 u,
          const std::vector<GridStructure<1>> &grid_struct,
          const std::vector<SolutionSnapshot<1>> &phi_history) const;
 

nufi/parameters.h

@@ -23,7 +23,7 @@ constexpr unsigned int NV = 128;
 constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV;
 
 // deal.ii options
-constexpr unsigned int GLOBAL_REFINEMENT = 8;
+constexpr unsigned int GLOBAL_REFINEMENT = 6;
 constexpr unsigned int FE_DEGREE = 3;
 constexpr unsigned int CONVERGENCE_ITERATIONS = 5000;
 constexpr double CONVERGENCE_LIMIT = 1e-8;
@@ -39,10 +39,10 @@ constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
 // NUFI options
 constexpr double DT = 1. / 10.;
 constexpr unsigned int TMAX = 100;
-constexpr unsigned int REFINE_FREQUENCY = 3;
+constexpr unsigned int REFINE_FREQUENCY = 30;
 
 // Plotting options
-constexpr int PLOT_FREQUENCY = 1;
+constexpr int PLOT_FREQUENCY = 10;
 constexpr size_t PLOT_NX = CALC_NX;
 constexpr double PLOT_DX = LX / PLOT_NX;
 const std::string PLOT_DIR = "results/";

nufi/poisson_problem.h

@@ -103,43 +103,27 @@ public:
 
   void save_grid_to_file(std::string &filename) const;
 
-  Triangulation<dim> triangulation;
-  DoFHandler<dim> dof_handler;
-  CellLocator<dim> cell_locator;
-
-  Vector<double> solution; // phi
-  SparsityPattern sparsity_pattern;
-  SparseMatrix<double> system_matrix;
-  Vector<double> system_rhs;
-
 private:
   void create_mesh();
   void setup_system();
   void assemble_system();
   void solve(size_t it);
 
-  // Triangulation<dim> triangulation;
-  FE_Q<dim> fe;
-  // DoFHandler<dim> dof_handler;
-
-  AffineConstraints<double> constraints;
-
-  // SparsityPattern sparsity_pattern;
-  // SparseMatrix<double> system_matrix;
-
-  // Vector<double> solution; // phi
-  // Vector<double> system_rhs;
-
   std::function<std::vector<double>(const std::vector<Point<dim>> &)>
       rhs_function;
-  // Vector<double> rhs;
 
   MappingQ<dim> mapping;
+  FE_Q<dim> fe;
+  AffineConstraints<double> constraints;
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler;
+  CellLocator<dim> cell_locator;
+  Vector<double> solution; // phi
+  SparsityPattern sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+  Vector<double> system_rhs;
 
   const bool PRINT_GAUGE_DOF_POSITION = true;
-
-  // mutable std::vector<double> local_solution_buffer;
-  // mutable std::unique_ptr<FEPointEvaluation<dim, dim>> evaluator;
 };
 
 //====//====//
@@ -164,10 +148,10 @@ void PoissonProblem<dim>::set_rhs_function(
 
 template <int dim>
 PoissonProblem<dim>::PoissonProblem(unsigned int degree)
-    : // triangulation(Triangulation<dim>::limit_level_difference_at_vertices),
-      triangulation(), dof_handler(triangulation), fe(degree), mapping(degree) {
+    : mapping(degree), fe(degree), triangulation(), dof_handler(triangulation) {
 }
 
+// Uses FEPointEvaluation
 template <int dim>
 std::vector<double>
 PoissonProblem<dim>::sample_electric_field(double x_min, double x_max,
@@ -180,32 +164,27 @@ PoissonProblem<dim>::sample_electric_field(double x_min, double x_max,
     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];
   }
 
@@ -471,6 +450,8 @@ template <int dim> void PoissonProblem<dim>::coarse_and_refine_grid(size_t it) {
 
 template <int dim> void PoissonProblem<dim>::solve(size_t it) {
 
+  std::cout << "Calling PoissonProblem::solve for time-step " << it << "\n";
+
   SolverControl solver_control(Parameters::CONVERGENCE_ITERATIONS,
                                Parameters::CONVERGENCE_LIMIT *
                                    system_rhs.l2_norm());
@@ -479,26 +460,26 @@ template <int dim> void PoissonProblem<dim>::solve(size_t it) {
   solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
   constraints.distribute(solution);
 
-  std::ofstream out("results/phi_after_solve_" + std::to_string(it) + ".dat");
-  std::vector<std::pair<double, double>> data;
-
-  const auto support =
-      DoFTools::map_dofs_to_support_points(mapping, dof_handler);
-
-  for (const auto &[dof, p] : support) {
-    data.emplace_back(p[0], solution[dof]);
-  }
-
-  std::sort(data.begin(), data.end());
-
-  for (const auto &[x, value] : data) {
-    out << x << " " << value << "\n";
-  }
-
-  std::vector<double> E_x =
-      sample_electric_field(Parameters::X_DOMAIN_LEFT,
-                            Parameters::X_DOMAIN_RIGHT, Parameters::PLOT_NX);
-  save_space_vector(E_x, "E_x_after_solve", it);
+  // std::ofstream out("results/phi_after_solve_" + std::to_string(it) +
+  // ".dat"); std::vector<std::pair<double, double>> data;
+  //
+  // const auto support =
+  //     DoFTools::map_dofs_to_support_points(mapping, dof_handler);
+  //
+  // for (const auto &[dof, p] : support) {
+  //   data.emplace_back(p[0], solution[dof]);
+  // }
+  //
+  // std::sort(data.begin(), data.end());
+  //
+  // for (const auto &[x, value] : data) {
+  //   out << x << " " << value << "\n";
+  // }
+  //
+  // std::vector<double> E_x =
+  //     sample_electric_field(Parameters::X_DOMAIN_LEFT,
+  //                           Parameters::X_DOMAIN_RIGHT, Parameters::PLOT_NX);
+  // save_space_vector(E_x, "E_x_after_solve", it);
 }
 
 template <int dim> void PoissonProblem<dim>::initialize() {
@@ -512,13 +493,10 @@ void PoissonProblem<dim>::solve_step(
   if (refining) {
     coarse_and_refine_grid(it);
     setup_system();
+    update_grid_versions(grid_versions, *this);
   }
-
   assemble_system();
   solve(it);
-
-  if (refining)
-    update_grid_versions(grid_versions, *this);
 }
 
 // NuFI doesnt use this, kept only for testing PoissonProblem

nufi/save_results.h

@@ -19,14 +19,9 @@ void save_rho(const NuFISolver &solver, unsigned int n,
               std::vector<SolutionSnapshot<1>> &phi_history,
               unsigned int Nx_out, const std::string &filename);
 
-void save_Efield(unsigned int n, GridStructure<1> &grid_struct,
+void save_Efield(unsigned int it, std::vector<GridStructure<1>> &grid_versions,
                  std::vector<SolutionSnapshot<1>> &phi_history,
-                 unsigned int Nx_out, const std::string &filename);
-
-void save_Efield_new(unsigned int it,
-                     std::vector<GridStructure<1>> &grid_versions,
-                     std::vector<SolutionSnapshot<1>> &phi_history,
-                     unsigned int Nx_out = Parameters::PLOT_NX);
+                 unsigned int Nx_out = Parameters::PLOT_NX);
 
 void save_space_vector(const std::vector<double> &vals,
                        const std::string &filename, size_t it);