Merge pull request #154 from AetherModel/run_in_1d

Code can run in 1D (alt)

Author: aaronjridley

include/aether.h

@@ -4,6 +4,8 @@
 #ifndef INCLUDE_AETHER_H_
 #define INCLUDE_AETHER_H_
 
+const float aether_version = 0.6;
+
 /// 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>

include/grid.h

@@ -172,6 +172,14 @@ class Grid {
 
   int64_t get_nGCs();
 
+  bool get_HasXdim();
+  bool get_HasYdim();
+  bool get_HasZdim();
+  bool get_Is0D();
+  bool get_Is1Dx();
+  bool get_Is1Dy();
+  bool get_Is1Dz();
+
   void fill_grid(Planets planet);
   void correct_xy_grid(Planets planet);
   void calc_sza(Planets planet, Times time);
@@ -311,6 +319,19 @@ class Grid {
   int64_t nY, nLats;
   int64_t nZ, nAlts;
 
+  // These logicals define the dimensionality of the grid:
+  bool Is0D = false;
+  bool Is1Dx = false;
+  bool Is1Dy = false;
+  bool Is1Dz = false;
+  bool Is2Dxy = false;
+  bool Is2Dxz = false;
+  bool Is2Dyz = false;
+  bool Is3D = false;
+  bool HasXdim = true;
+  bool HasYdim = true;
+  bool HasZdim = true;
+
   int nGCs; // number of ghostcells
 
   // interpolation members

include/output.h

@@ -119,6 +119,12 @@ class OutputContainer {
    **/
   void set_version(float in_version);
 
+  /**********************************************************************
+     \brief set the number of ghost cells in Aether
+     \param in_nGCs the number of ghostcells in the grid
+   **/
+  void set_nGhostCells(int in_nGCs);
+
   /**********************************************************************
      \brief write a file with the information in the container
    **/
@@ -191,9 +197,13 @@ class OutputContainer {
 
   /// The time of the data
   std::vector<int> itime;
+
   /// The version of the code / data / whatever:
   float version;
 
+  /// The number of Ghost Cells:
+  int nGCs;
+
   /// The frequency of the output for this particular container:
   float dt_output;
 

share/run/UA/inputs/defaults.json

@@ -69,8 +69,8 @@
 	"LatRange" : [-90.0, 90.0],
         "nLatsPerBlock" : 18,
 	"LonRange" : [0.0, 360.0],
-        "nLonsPerBlock" : 36,
-        "nAlts" : 50,
+        "nLonsPerBlock" : 20,
+        "nAlts" : 40,
 	"MinAlt" : 100.0,
 	"dAltkm" : 5.0,
         "dAltScale" : 0.25,
@@ -82,11 +82,15 @@
 	"LatRange" : [-90.0, 90.0],
         "nLatsPerBlock" : 18,
 	"LonRange" : [0.0, 360.0],
-        "nLonsPerBlock" : 36,
-        "nAlts" : 200,
+        "nLonsPerBlock" : 22,
+        "nAlts" : 50,
 	"MinAlt" : 80.0,
         "MinApex" : 120.0,
-        "MaxAlt" : 5000.0},
+        "MaxAlt" : 5000.0,
+	"dAltkm" : 5.0,
+        "dAltScale" : 0.25,
+	"IsUniformAlt" : true,
+        "AltFile" : ""},
 
     "Oblate" : {
 	"isOblate" : false,

src/containers_io.cpp

@@ -162,6 +162,14 @@ void OutputContainer::set_version(float in_version) {
   version = in_version;
 }
 
+// -----------------------------------------------------------------------------
+// Set the number of ghostcells.
+// -----------------------------------------------------------------------------
+
+void OutputContainer::set_nGhostCells(int in_nGCs) {
+  nGCs = in_nGCs;
+}
+
 // -----------------------------------------------------------------------------
 // Clears the elements vector within the output container
 // -----------------------------------------------------------------------------
@@ -176,6 +184,7 @@ void OutputContainer::clear_variables() {
 
 OutputContainer::OutputContainer() {
   // Set default output type to netCDF
+  nGCs = 0;
 #ifdef NETCDF
   output_type = netcdf_type;
 #else
@@ -236,6 +245,7 @@ void OutputContainer::display() {
   std::cout << "  nX : " << elements[0].value.n_rows << "\n";
   std::cout << "  nY : " << elements[0].value.n_cols << "\n";
   std::cout << "  nZ : " << elements[0].value.n_slices << "\n";
+  std::cout << "  nGCs : " << nGCs << "\n";
   int64_t nVars = elements.size();
   std::cout << "  Number of Variables : " << nVars << "\n";
 
@@ -281,6 +291,7 @@ bool OutputContainer::write_container_header() {
     {"nX", nX},
     {"nY", nY},
     {"nZ", nZ},
+    {"nGCs", nGCs},
     {"nLons", nX},
     {"nLats", nY},
     {"nAlts", nZ},
@@ -327,6 +338,7 @@ bool OutputContainer::read_container_binary() {
   int64_t iY, nY = header["nY"];
   int64_t iZ, nZ = header["nZ"];
   int64_t iTotalSize = nX * nY * nZ;
+  nGCs = header["nGCs"];
 
   float *variable_array = new float[iTotalSize];
   arma_cube value_scgc;

src/exchange_messages.cpp

@@ -896,6 +896,10 @@ bool exchange_one_var(Grid &grid,
                       arma_cube &var_to_pass,
                       bool doReverseSignAcrossPole) {
 
+  // This function is only needed if we do interpolation, which only happens in
+  // the horizontal directions
+  if (!grid.get_HasXdim() & !grid.get_HasYdim()) return true;
+
   std::string function = "exchange_one_var";
   static int iFunction = -1;
   report.enter(function, iFunction);
@@ -1192,6 +1196,10 @@ arma_cube interpolate_ghostcells(arma_cube varIn, Grid &grid) {
 
 bool find_ghostcell_interpolation_coefs(Grid &grid) {
 
+  // This function is only needed if we do interpolation, which only happens in
+  // the horizontal directions
+  if (!grid.get_HasXdim() & !grid.get_HasYdim()) return true;
+
   std::string function = "find_ghostcell_interpolation_coefs";
   static int iFunction = -1;
   report.enter(function, iFunction);

src/grid.cpp

@@ -23,6 +23,28 @@ Grid::Grid(std::string gridtype) {
   nZ = grid_input.nZ + nGCs * 2;
   nAlts = nZ;
 
+  // No set all of the logicals to make the flow a bit easier:
+
+  if (grid_input.nX == 1 &
+      grid_input.nY == 1 &
+      grid_input.nZ == 1)
+    Is0D = true;
+  else {
+    if (grid_input.nY == 1 & grid_input.nZ == 1) Is1Dx = true;
+    if (grid_input.nX == 1 & grid_input.nZ == 1) Is1Dy = true;
+    if (grid_input.nX == 1 & grid_input.nY == 1) Is1Dz = true;
+    if (!Is1Dx & !Is1Dy & !Is1Dz) {
+      if (grid_input.nX == 1) Is2Dyz = true;
+      if (grid_input.nY == 1) Is2Dxz = true;
+      if (grid_input.nZ == 1) Is2Dxy = true;
+      if (!Is2Dyz & !Is2Dxz & !Is2Dxy) Is3D = true;
+    }
+  }
+
+  if (grid_input.nX == 1) HasXdim = false;
+  if (grid_input.nY == 1) HasYdim = false;
+  if (grid_input.nZ == 1) HasZdim = false;
+
   if (mklower(grid_input.shape) == "sphere") 
     iGridShape_ = iSphere_;
   if (mklower(grid_input.shape) == "cubesphere") 
@@ -192,8 +214,9 @@ bool Grid::write_restart(std::string dir) {
     try {
       OutputContainer RestartContainer;
       RestartContainer.set_directory(dir);
-      RestartContainer.set_version(0.1);
+      RestartContainer.set_version(aether_version);
       RestartContainer.set_time(0.0);
+      RestartContainer.set_nGhostCells(nGCs);
 
       // Output Cell Centers
       RestartContainer.set_filename("grid_" + cGrid);
@@ -280,7 +303,7 @@ bool Grid::read_restart(std::string dir) {
   try {
     OutputContainer RestartContainer;
     RestartContainer.set_directory(dir);
-    RestartContainer.set_version(0.1);
+    RestartContainer.set_version(aether_version);
     // Cell Centers:
     RestartContainer.set_filename("grid_" + cGrid);
     RestartContainer.read();
@@ -402,6 +425,38 @@ int64_t Grid::get_nPointsInGrid() {
   return nPoints;
 }
 
+// --------------------------------------------------------------------------
+// Get some grid definition things
+// --------------------------------------------------------------------------
+
+bool Grid::get_HasXdim() {
+  return HasXdim;
+}
+
+bool Grid::get_HasYdim() {
+  return HasYdim;
+}
+
+bool Grid::get_HasZdim() {
+  return HasZdim;
+}
+
+bool Grid::get_Is0D() {
+  return Is0D;
+}
+
+bool Grid::get_Is1Dx() {
+  return Is1Dx;
+}
+
+bool Grid::get_Is1Dy() {
+  return Is1Dy;
+}
+
+bool Grid::get_Is1Dz() {
+  return Is1Dz;
+}
+
 // --------------------------------------------------------------------------
 // Get a bunch of sizes within the grid
 // --------------------------------------------------------------------------

src/init_geo_grid.cpp

@@ -650,6 +650,10 @@ void Grid::create_sphere_grid(Quadtree quadtree) {
   precision_t lon0 = lower_left_norm(0) * cPI;
   arma_vec lon1d(nLons);
 
+  // if we are not doing anything in the lon direction, then set dlon to
+  // something reasonable:
+  if (!HasXdim) dlon = 1.0 * cDtoR;
+
   // Longitudes:
   // - Make a 1d vector
   // - copy it into the 3d cube
@@ -665,6 +669,10 @@ void Grid::create_sphere_grid(Quadtree quadtree) {
   precision_t lat0 = lower_left_norm(1) * cPI;
   arma_vec lat1d(nLats);
 
+  // if we are not doing anything in the lat direction, then set dlat to
+  // something reasonable:
+  if (!HasYdim) dlat = 1.0 * cDtoR;
+
   // Latitudes:
   // - Make a 1d vector
   // - copy it into the 3d cube
@@ -947,11 +955,11 @@ bool Grid::init_geo_grid(Quadtree quadtree,
 
   if (iGridShape_ == iCubesphere_) {
     report.print(0, "Creating Cubesphere Grid");
-    create_cubesphere_connection(quadtree);
+    if (!Is0D & !Is1Dz) create_cubesphere_connection(quadtree);
     IsCubeSphereGrid = true;
   } else {
     report.print(0, "Creating Spherical Grid");
-    create_sphere_connection(quadtree);
+    if (!Is0D & !Is1Dz) create_sphere_connection(quadtree);
     IsCubeSphereGrid = false;
   }
 

src/ions.cpp

@@ -303,7 +303,7 @@ bool Ions::restart_file(std::string dir, bool DoRead) {
     if (DoRead)
       RestartContainer.read();
     else {
-      RestartContainer.set_version(0.1);
+      RestartContainer.set_version(aether_version);
       RestartContainer.set_time(0.0);
     }
 

src/neutrals.cpp

@@ -359,7 +359,7 @@ bool Neutrals::restart_file(std::string dir, bool DoRead) {
     if (DoRead)
       RestartContainer.read();
     else {
-      RestartContainer.set_version(0.1);
+      RestartContainer.set_version(aether_version);
       RestartContainer.set_time(0.0);
     }