grid refinement on, eval needs update to use new grid, solver not converging on time

nufi/parameters.h

@@ -3,6 +3,7 @@
 
 #include <cmath>
 #include <cstdlib>
+#include <string>
 
 namespace Parameters {
 constexpr unsigned int DIMENSION = 1;
@@ -24,7 +25,7 @@ constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV;
 // deal.ii options
 constexpr unsigned int GLOBAL_REFINEMENT = 6;
 constexpr unsigned int FE_DEGREE = 2;
-constexpr unsigned int CONVERGENCE_ITERATIONS = 5000;
+constexpr unsigned int CONVERGENCE_ITERATIONS = 10000;
 constexpr double CONVERGENCE_LIMIT = 1e-8;
 
 constexpr double EPS = 0.01;
@@ -32,13 +33,15 @@ constexpr double WAVE_NR = 0.5;
 constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
 
 // NUFI options
-constexpr double DT = 1. / 8.;
+constexpr double DT = 1. / 10.;
 constexpr unsigned int TMAX = 100;
+constexpr unsigned int REFINE_FREQUENCY = 5;
 
 // Plotting options
-constexpr int PLOT_FREQUENCY = 4;
+constexpr int PLOT_FREQUENCY = 5;
 constexpr size_t PLOT_NX = CALC_NX;
 constexpr double PLOT_DX = LX / PLOT_NX;
+const std::string PLOT_DIR = "results/";
 } // namespace Parameters
 
 #endif

nufi/poisson_problem.h

@@ -66,7 +66,8 @@ public:
 
   void initialize();
   void solve_step();
-  void coarse_and_refine_grid();
+  void coarse_and_refine_grid(size_t it,
+                              std::vector<Vector<double>> &solution_history);
   void run();
 
   void set_rhs_function(std::function<double(const Point<dim> &)> f);
@@ -84,7 +85,10 @@ public:
   eval_vector_grad(const Vector<double> &solution,
                    const std::vector<Point<dim>> &points) const;
 
-  void save_grid_to_file(const std::string &filename) const;
+  void save_grid_to_file(std::string &filename) const;
+
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler;
 
 private:
   void create_mesh();
@@ -92,9 +96,9 @@ private:
   void assemble_system();
   void solve();
 
-  Triangulation<dim> triangulation;
+  // Triangulation<dim> triangulation;
   FE_Q<dim> fe;
-  DoFHandler<dim> dof_handler;
+  // DoFHandler<dim> dof_handler;
 
   AffineConstraints<double> constraints;
 
@@ -250,10 +254,18 @@ double eval_point_value(const Mapping<dim> &mapping,
 }
 
 template <int dim>
-void PoissonProblem<dim>::save_grid_to_file(const std::string &filename) const {
+void PoissonProblem<dim>::save_grid_to_file(std::string &filename) const {
-  std::ofstream out(filename);
   GridOut grid_out;
-  grid_out.write_svg(triangulation, out);
+
+  if (dim >= 2) {
+    filename += ".svg";
+    std::ofstream out(filename);
+    grid_out.write_svg(triangulation, out);
+  } else if (dim == 1) {
+    filename += ".gnuplot";
+    std::ofstream out(filename);
+    grid_out.write_gnuplot(triangulation, out);
+  }
 
   std::cout << "Grid written to " << filename << "\n";
 }
@@ -371,12 +383,8 @@ template <int dim> void PoissonProblem<dim>::assemble_system() {
 
 template <int dim>
 void PoissonProblem<dim>::coarse_and_refine_grid(
-    std::vector<Vector<double>> &solution_history) {
+    size_t it, std::vector<Vector<double>> &solution_history) {
-  // Add refinement and coasring algorithm here
+  std::cout << "Refinement Started" << "\n";
-
-  //======//======//
-  // CHECK step-15.
-  //======//======//
   Vector<float> error_per_cell(triangulation.n_active_cells());
 
   KellyErrorEstimator<dim>::estimate(
@@ -384,15 +392,10 @@ void PoissonProblem<dim>::coarse_and_refine_grid(
       std::map<types::boundary_id, const Function<dim> *>(), solution,
       error_per_cell);
 
-  GridRefinement::refine_and_coarsen_fixed_number(
+  GridRefinement::refine_and_coarsen_fixed_number(triangulation, error_per_cell,
-      triangulation, error_per_cell, 0.3,
+                                                  0.3, 0.03);
-      0.03); // These numbers are to be reviewd and changed
 
   triangulation.prepare_coarsening_and_refinement();
-
-  // old solutions transfer
-  const size_t N_sols = solution_history.size();
-
   SolutionTransfer<dim, Vector<double>> transfer(dof_handler);
   transfer.prepare_for_coarsening_and_refinement(solution_history);
 
@@ -400,14 +403,27 @@ void PoissonProblem<dim>::coarse_and_refine_grid(
 
   setup_system();
 
-  transfer.interpolate(solution_history);
+  std::vector<Vector<double>> new_solution_history(solution_history.size());
+
+  for (auto &vec : new_solution_history)
+    vec.reinit(dof_handler.n_dofs());
+
+  transfer.interpolate(solution_history, new_solution_history);
+
+  solution_history.swap(new_solution_history);
 
   solution = solution_history.back();
 
   constraints.distribute(solution);
 
-  // rebuild cells with the grid
   cell_locator.rebuild(dof_handler, triangulation);
+
+  std::cout << "Refinement Finished" << "\n";
+
+  std::string grid_file_name =
+      Parameters::PLOT_DIR + "grid_" + std::to_string(it);
+
+  save_grid_to_file(grid_file_name);
 }
 
 template <int dim> void PoissonProblem<dim>::solve() {

src/nufi_solver.cc

@@ -173,7 +173,14 @@ void NuFISolver::run() {
   std::ofstream time_file("results/simulation_time.dat");
 
   double total_time = 0;
-  time_file << "it step_time total_time" << "\n";
+
+  time_file << "it "
+            << "step_time "
+            << "total_time "
+            << "compute_time "
+            << "refine_time "
+            << "plot_time"
+            << "\n";
 
   const double x_min = Parameters::X_DOMAIN_LEFT;
   double dx = Parameters::CALC_DX;
@@ -182,11 +189,30 @@ void NuFISolver::run() {
     stopwatch<double> timer;
 
     double time_elapsed_before = timer.elapsed();
+    double compute_time = 0.0;
+    double refine_time = 0.0;
+    double plot_time = 0.0;
 
     std::cout << "Timestep " << it << " / " << Nt
               << " (simulation time = " << it * Parameters::DT << ")"
               << std::endl;
 
+    // START: diagnostics
+    std::cout << "cells = " << poisson.triangulation.n_active_cells()
+              << " dofs = " << poisson.dof_handler.n_dofs() << std::endl;
+    double min_h = 1e100;
+    double max_h = 0;
+
+    for (auto cell : poisson.triangulation.active_cell_iterators()) {
+      min_h = std::min(min_h, cell->diameter());
+      max_h = std::max(max_h, cell->diameter());
+    }
+
+    std::cout << "h ratio = " << max_h / min_h << std::endl;
+
+    // END: diagnostics
+
+    double compute_start = timer.elapsed();
     // compute rho
     std::vector<double> x_eval = make_x_eval(Nx);
     std::vector<double> tmp_rho =
@@ -207,19 +233,20 @@ void NuFISolver::run() {
     poisson.solve_step();
     phi_history.push_back(poisson.get_solution());
 
-    // std::vector<double> sampled_potential =
+    compute_time = timer.elapsed() - compute_start;
-    //     poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of
+
-    //     FE
+    if (it % Parameters::REFINE_FREQUENCY == 0 && it != 0) {
+      double refine_start = timer.elapsed();
+      poisson.coarse_and_refine_grid(it, phi_history);
+      refine_time = timer.elapsed() - refine_start;
+    }
 
     double timer_elapsed = timer.elapsed();
     double step_time = timer_elapsed - time_elapsed_before;
-    total_time += timer_elapsed;
-
-    time_file << it << " " << step_time << " " << total_time << "\n";
-    time_file.flush();
 
     std::cout << "step made in " << step_time << " seconds\n\n";
     if (it % Parameters::PLOT_FREQUENCY == 0) {
+      double plot_start = timer.elapsed();
       std::cout << "Saving results...   ";
       save_f(*this, it, poisson, phi_history, Parameters::PLOT_NX,
              Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
@@ -241,7 +268,13 @@ void NuFISolver::run() {
       int_E_squared.push_back(int_val);
       save_space_vector(int_E_squared, "electricint", it);
       std::cout << "Time since start = " << total_time << "\n\n";
+
+      plot_time = timer.elapsed() - plot_start;
     }
+
+    time_file << it << " " << step_time << " " << total_time << " "
+              << compute_time << " " << refine_time << " " << plot_time << "\n";
+    time_file.flush();
   }
 
   std::cout << "NuFI simulation finished in " << total_time << " seconds.\n";