Test vertical advect

edu/examples/Advection/advect.cpp

@@ -1,7 +1,31 @@
 
 // g++ -I/usr/local/include -I/Users/ridley/Software/Json/json/include main.cpp
+// g++ -I/usr/local/include -o advect1d advect.cpp
+
+/// The armadillo library is to allow the use of 3d cubes and other
+/// array types, with array math built in. This eliminates loops!
+#include <armadillo>
+
+/// This is used for timing and the random seed generator:
+#include <chrono>
+
+// Types
+// Precision compile-time aliasing
+#ifdef AETHER_USE_PRECISION_DOUBLE
+/// Precision type chosen to be `double` through `AETHER_USE_PRECISION_DOUBLE`
+using precision_t = double;
+#else
+/// Precision type compile-time default to float.
+using precision_t = float;
+#endif
+
+/// Armadillo type vector (single column) with compile-time precision.
+using arma_vec = arma::Col<precision_t>;
+/// Armadillo type matrix (two dimension) with compile-time precision.
+using arma_mat = arma::Mat<precision_t>;
+/// Armadillo type cube (three dimension) with compile-time precision.
+using arma_cube = arma::Cube<precision_t>;
 
-#include "../../../include/aether.h"
 #include <fstream>
 
 // ---------------------------------------------------------
@@ -17,10 +41,40 @@ arma_vec init_grid(int64_t nPts, int64_t nGCs) {
   for (int64_t i = -nGCs; i < nPts + nGCs; i++) {
     x(i + nGCs) = i * dx;
   }
+  precision_t maxX = x(nPts + nGCs - 1);
+  x = 100.0 * x / maxX;
 
   return x;
 }
 
+// ---------------------------------------------------------
+// grid stretched creation
+// ---------------------------------------------------------
+
+arma_vec init_stretched_grid(int64_t nPts, int64_t nGCs) {
+
+  precision_t dx = 1.0;
+  arma_vec x(nPts + nGCs * 2);
+
+  precision_t factor = 1.0;
+  precision_t i2pi = 2.0 * 3.1415927 / (nPts-1);
+
+  x(nGCs) = 0.0;
+
+  for (int64_t i = 1; i < nPts + nGCs; i++) {
+    x(i + nGCs) = x(i - 1 + nGCs) + dx  + factor * (1 + cos(i * i2pi));
+    std::cout << "i : " << i << " " << cos(i * i2pi) << "\n";
+  }
+  for (int64_t i = -1; i >= -nGCs; i--) {
+    x(i + nGCs) = x(i + 1 + nGCs) - dx  - factor * (1 + cos(i * i2pi));
+    std::cout << "i : " << i << " " << cos(i * i2pi) << "\n";
+  }
+  precision_t maxX = x(nPts + nGCs - 1);
+  x = 100.0 * x / maxX;
+
+  return x;
+}
+
 // ---------------------------------------------------------
 // bin edges
 // ---------------------------------------------------------
@@ -80,6 +134,7 @@ arma_vec init_vel(int64_t nPts) {
   arma_vec vel(nPts);
   // all cells positive to right:
   vel.ones();
+  vel = -1.0 * vel;
   return vel;
 }
 
@@ -182,6 +237,43 @@ arma_vec calc_grad(arma_vec values,
   return gradients;
 }
 
+// ---------------------------------------------------------
+// Limiter on values
+//   projected is assumed to be on the edge between the
+//   i-1 and i cell (i-1/2)
+//   limited is returned at edges
+// ---------------------------------------------------------
+
+arma_vec limiter_value(arma_vec projected,
+		       arma_vec values,
+		       int64_t nPts,
+		       int64_t nGCs) {
+
+  int64_t iStart = 0;
+  int64_t iEnd = nPts + 2 * nGCs;
+
+  arma_vec limited = projected;
+
+  precision_t mini, maxi;
+
+  for (int64_t i = iStart + 1; i < iEnd - 1; i++) {
+
+    mini = values(i-1);
+    if (values(i) < mini)
+      mini = values(i);
+    maxi = values(i-1);
+    if (values(i) > maxi)
+      maxi = values(i);
+
+    if (limited(i) < mini)
+      limited(i) = mini;
+    if (limited(i) > maxi)
+      limited(i) = maxi;
+
+  }
+  return limited;
+}
+
 // ---------------------------------------------------------
 // Project gradients + values to the right face, from the left
 //   returned values are on the i - 1/2 edges
@@ -209,6 +301,42 @@ arma_vec project_from_left(arma_vec values,
   return projected;
 }
 
+// ---------------------------------------------------------
+// Project gradients + values to the right face, from the left
+//   returned values are on the i - 1/2 edges
+//     (between i-1 and i cell center)
+// ---------------------------------------------------------
+
+arma_vec project_from_left_new(arma_vec values,
+			       arma_vec x_centers,
+			       arma_vec x_edges,
+			       int64_t nPts,
+			       int64_t nGCs) {
+  int64_t iStart = 1;
+  int64_t iEnd = nPts + 2 * nGCs - 1;
+
+  // Define at edges:
+  arma_vec projected(nPts + 2 * nGCs + 1);
+  projected.zeros();
+
+  precision_t dxei, dxci, dxcip1, r;
+
+  // no gradient in the 0 or iEnd cells
+  for (int64_t i = iStart; i < iEnd; i++) {
+    dxei = x_edges(i + 1) - x_edges(i);
+    dxci = x_centers(i) - x_centers(i - 1);
+    dxcip1 = x_centers(i + 1) - x_centers(i);
+    r = dxcip1 / dxci;
+    projected(i + 1) = values(i) + 
+      0.5 * dxei * (values(i) - values(i - 1)) / dxci + 
+      0.125 * dxei * dxei * (values(i + 1) + r * values(i - 1) - (1 + r) * values(i)) / (dxci * dxcip1);
+  }
+
+  projected = limiter_value(projected, values, nPts, nGCs);
+
+  return projected;
+}
+
 
 // ---------------------------------------------------------
 // Project gradients + values to the left face, from the right
@@ -238,40 +366,39 @@ arma_vec project_from_right(arma_vec values,
 }
 
 // ---------------------------------------------------------
-// Limiter on values
-//   projected is assumed to be on the edge between the
-//   i-1 and i cell (i-1/2)
-//   limited is returned at edges
+// Project gradients + values to the left face, from the right
+//   returned values are on the i - 1 edges
+//     (between i-1 and i cell center)
 // ---------------------------------------------------------
 
-arma_vec limiter_value(arma_vec projected,
-		       arma_vec values,
-		       int64_t nPts,
-		       int64_t nGCs) {
-
-  int64_t iStart = 0;
-  int64_t iEnd = nPts + 2 * nGCs;
-
-  arma_vec limited = projected;
+arma_vec project_from_right_new(arma_vec values,
+				arma_vec x_centers,
+				arma_vec x_edges,
+				int64_t nPts,
+				int64_t nGCs) {
+  int64_t iStart = 1;
+  int64_t iEnd = nPts + 2 * nGCs - 1;
 
-  precision_t mini, maxi;
+  // Define at edges:
+  arma_vec projected(nPts + 2 * nGCs + 1);
+  precision_t dxei, dxci, dxcip1, r;
 
-  for (int64_t i = iStart + 1; i < iEnd - 1; i++) {
-
-    mini = values(i-1);
-    if (values(i) < mini)
-      mini = values(i);
-    maxi = values(i-1);
-    if (values(i) > maxi)
-      maxi = values(i);
+  projected.zeros();
 
-    if (limited(i) < mini)
-      limited(i) = mini;
-    if (limited(i) > maxi)
-      limited(i) = maxi;
-
+  // no gradient in the 0 or iEnd cells
+  for (int64_t i = iStart; i < iEnd; i++) {
+    dxei = x_edges(i + 1) - x_edges(i);
+    dxci = x_centers(i) - x_centers(i - 1);
+    dxcip1 = x_centers(i + 1) - x_centers(i);
+    r = dxcip1 / dxci;
+    projected(i) = values(i) - 
+      0.5 * dxei * (values(i + 1) - values(i)) / dxcip1 +
+      0.125 * dxei * dxei * (values(i + 1) + r * values(i - 1) - (1 + r) * values(i)) / (dxci * dxcip1);
   }
-  return limited;
+
+  projected = limiter_value(projected, values, nPts, nGCs);
+
+  return projected;
 }
 
 // ---------------------------------------------------------
@@ -363,19 +490,21 @@ void output(arma_vec values,
 
 int main() {
 
-  precision_t dt = 0.01;
-
-  int64_t nSteps = 3;
-  int64_t iStep;
-
-  int64_t nPts = 60;
+  int64_t nPts = 200;
   int64_t nGCs = 2;
   int64_t nPtsTotal = nGCs + nPts + nGCs;
 
   arma_vec x = init_grid(nPts, nGCs);
+  //arma_vec x = init_stretched_grid(nPts, nGCs);
   arma_vec edges = calc_bin_edges(x);
   arma_vec widths = calc_bin_widths(edges);
 
+  precision_t dt = 0.1 * x(nPts + nGCs - 1) / nPts;
+  precision_t time = 0.0;
+
+  int64_t nSteps = x(nPts + nGCs - 1) / dt;
+  int64_t iStep;
+
   arma_vec grad_rho;
   arma_vec rhoL;
   arma_vec rhoR;
@@ -392,19 +521,29 @@ int main() {
   exchange(vel, nPts, nGCs);
 
   output(rho, "rho.txt", false, nPts, nGCs);
+  output(x, "x.txt", false, nPts, nGCs);
 
   for (iStep = 0; iStep < nSteps; iStep++) {
+
+    std::cout << "iStep = " << iStep << "; time =  " << time << "\n";
+    time = time + dt;
 
     grad_rho = calc_grad(rho, x, nPts, nGCs);
 
     // Right side of edge from left
-    rhoR = project_from_left(rho, grad_rho,
+//rhoR = project_from_left(rho, grad_rho,
+//			     x, edges,
+//			     nPts, nGCs);
+    rhoR = project_from_left_new(rho,
 			     x, edges,
 			     nPts, nGCs);
     //rhoR = limiter_value(rhoR, rho, nPts, nGCs);
 
     // Left side of edge from left
-    rhoL = project_from_right(rho, grad_rho,
+//    rhoL = project_from_right(rho, grad_rho,
+//			      x, edges,
+//			      nPts, nGCs);
+    rhoL = project_from_right_new(rho,
 			      x, edges,
 			      nPts, nGCs);
     //rhoL = limiter_value(rhoL, rho, nPts, nGCs);
@@ -413,13 +552,19 @@ int main() {
     grad_vel = calc_grad(vel, x, nPts, nGCs);
 
     // Right side of edge from left
-    velR = project_from_left(vel, grad_vel,
+//    velR = project_from_left(vel, grad_vel,
+//			     x, edges,
+//			     nPts, nGCs);
+    velR = project_from_left_new(vel,
 			     x, edges,
 			     nPts, nGCs);
     //velR = limiter_value(velR, vel, nPts, nGCs);
 
     // Left side of edge from left
-    velL = project_from_right(vel, grad_vel,
+//    velL = project_from_right(vel, grad_vel,
+//			      x, edges,
+//			      nPts, nGCs);
+    velL = project_from_right_new(vel,
 			      x, edges,
 			      nPts, nGCs);
     //velL = limiter_value(velL, vel, nPts, nGCs);