- working solver

experiment.cpp

@@ -9,6 +9,16 @@ double index_to_x(std::size_t i, double xl, double dx)
   return xl + i * dx;
 }
 
+double flux(double u)
+{
+  return 2 * u;
+}
+
+double numflux(double ul, double ur)
+{
+  return flux(ul) + flux(ur) - 2 * (ur - ul);
+}
+
 int main(int argc, char** argv)
 {
   // Initialize the MPI environment
@@ -35,6 +45,7 @@ int main(int argc, char** argv)
   MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
 
   std::vector<double> data(ncell);
+  std::vector<double> tmp(ncell);
   std::vector<double> halo(nhalo);
 
   if (world_rank == 0) {
@@ -49,12 +60,38 @@ int main(int argc, char** argv)
     }
   }
 
-  auto const filename = "data_" + std::to_string(world_rank) + ".txt";
+  auto const filename = "data_" + std::to_string(world_rank) + ".0.txt";
   std::ofstream file(filename);
   for (std::size_t i = 0; i < ncell; ++i) {
     file << index_to_x(i, xl, dx) << " " << data[i] << "\n";
   }
 
+  auto dt = dx / 8;
+
+  double time = 0;
+  double Tfinal = 0.2;
+  std::size_t iter = 0;
+  std::size_t max_iter = 1000;
+
+  while (time < Tfinal and iter < max_iter) {
+    dt = std::min(dt, Tfinal - time);
+    for (std::size_t i = 1; i < ncell - 1; ++i) {
+      tmp[i] = data[i] + dt / dx * (numflux(data[i - 1], data[i]) - numflux(data[i], data[i + 1]));
+    }
+    for (std::size_t i = 1; i < ncell - 1; ++i) {
+      data[i] = tmp[i];
+    }
+    time += dt;
+    ++iter;
+    std::cout << "iter = " << iter << "\ttime = " << time << "\tdt = " << dt << "\n";
+  }
+
+  auto const final_filename = "data_" + std::to_string(world_rank) + ".final.txt";
+  std::ofstream final_file(final_filename);
+  for (std::size_t i = 0; i < ncell; ++i) {
+    final_file << index_to_x(i, xl, dx) << " " << data[i] << "\n";
+  }
+
   // Finalize the MPI environment.
   MPI_Finalize();
 }