Merge pull request #184 from AetherModel/messagepass_grids

This is for message passing between two grids within Aether.

Author: aaronjridley

include/calc_grid_derived.h

@@ -19,9 +19,14 @@ arma_vec calc_bin_widths(arma_vec centers);
 // ----------------------------------------------------------------------------
 // A helper function for mapping grids
 // ----------------------------------------------------------------------------
-bool grid_match(Grid gGrid,
-                Grid mGrid,
+bool grid_match(Grid &gGrid, 
+                Grid &mGrid, 
                 Quadtree gQuadtree,
                 Quadtree mQuadtree);
 
+bool get_data_from_other_grid(Grid &gGrid, 
+                              Grid &mGrid,
+                              arma_cube &gData,
+                              arma_cube &mData);
+
 #endif  // INCLUDE_CALC_GRID_DERIVED_H_

include/grid.h

@@ -45,17 +45,53 @@ struct cubesphere_chars {
 // Grid class
 // ----------------------------------------------------------------------------
 
-class Grid {
+struct interp_coef_t {
+  // The point is inside the cube of:
+  // [iRow, iRow+1], [iCol, iCol+1] [iAlt, iAlt+1]
+  uint64_t iRow;
+  uint64_t iCol;
+  uint64_t iAlt;
+  // The coefficients along row, column and altitude
+  precision_t rRow;
+  precision_t rCol;
+  precision_t rAlt;
+  // Whether the point is within this grid or not
+  bool in_grid;
+  // If this is set to true:
+  bool above_grid, below_grid;
+  // do interpolation in lat and lon, but extrapolate in altitude
+};
+
+struct grid_to_grid_t {
+  int64_t iProcTo;
+  int64_t nPts;
+  int64_t nPtsReceive;
+  std::vector<struct interp_coef_t> interpCoefs;
+  std::vector<precision_t *> valueToSend;
+  std::vector<precision_t *> valueToReceive;
+};
 
+class Grid {
  public:
   int iGridShape_ = -1;
-  // Armidillo Cube Versions:
+  // The index and coefficient used for interpolation
+  // Each point is processed by the function set_interpolation_coefs and stored
+  // in the form of this structure.
+  // If the point is out of the grid, in_grid = false and all other members are undefined
+
+  std::vector<struct grid_to_grid_t> gridToGridCoefs;
+  arma::Cube<int> gridToGridMap;
+
+  // Armadillo Cube Versions:
   // Cell Center Coordinates
   arma_cube geoLon_scgc, geoX_scgc;
   arma_cube geoLat_scgc, geoY_scgc;
   arma_cube geoAlt_scgc, geoZ_scgc;
   arma_cube geoLocalTime_scgc;
 
+  // This is an array for testing things:
+  arma_cube test_scgc;
+
   // Reference coordinate
   arma_cube refx_scgc, refy_scgc;
 
@@ -508,6 +544,21 @@ class Grid {
                                bool areLocsGeo = true,
                                bool areLocsIJK = true);
 
+  /**
+   * \brief Set the interpolation coefficients
+   * \param Lons The longitude of points
+   * \param Lats The latitude of points
+   * \param Alts The altitude of points
+   * \pre This instance is an geo grid
+   * \pre Lons, Lats and Alts have the same size
+   * \return list of interpolation coefficients
+   */
+
+  std::vector<struct interp_coef_t> get_interpolation_coefs(
+    const std::vector<precision_t> &Lons,
+    const std::vector<precision_t> &Lats,
+    const std::vector<precision_t> &Alts);
+
   /**
    * \brief Set the interpolation coefficients for the dipole grid
    * \param Lons The longitude of points
@@ -530,6 +581,8 @@ class Grid {
    *         an empty vector if the data is not the same size as the geo grid.
    */
   std::vector<precision_t> get_interpolation_values(const arma_cube &data) const;
+  std::vector<precision_t> get_interpolation_values(arma_cube data,
+                                                    std::vector<struct interp_coef_t> coefArray);
 
  private:
   bool IsGeoGrid;
@@ -608,22 +661,24 @@ class Grid {
     precision_t alt_max;
   };
 
+  // Return the index of the last element that has altitude smaller than or euqal to the input
+  uint64_t search_altitude(const precision_t alt_in) const;
   // The index and coefficient used for interpolation
   // Each point is processed by the function set_interpolation_coefs and stored
   // in the form of this structure.
   // If the point is out of the grid, in_grid = false and all other members are undefined
-  struct interp_coef_t {
-    // The point is inside the cube of [iRow, iRow+1], [iCol, iCol+1], [iAlt, iAlt+1]
-    uint64_t iRow;
-    uint64_t iCol;
-    uint64_t iAlt;
-    // The coefficients along row, column and altitude
-    precision_t rRow;
-    precision_t rCol;
-    precision_t rAlt;
-    // Whether the point is within this grid or not
-    bool in_grid;
-  };
+  //struct interp_coef_t {
+  //  // The point is inside the cube of [iRow, iRow+1], [iCol, iCol+1], [iAlt, iAlt+1]
+  //  uint64_t iRow;
+  //  uint64_t iCol;
+  //  uint64_t iAlt;
+  //  // The coefficients along row, column and altitude
+  //  precision_t rRow;
+  //  precision_t rCol;
+  //  precision_t rAlt;
+  //  // Whether the point is within this grid or not
+  //  bool in_grid;
+  //};
 
   // Calculate the range of a spherical grid
   void get_sphere_grid_range(struct sphere_range &sr) const;
@@ -633,19 +688,19 @@ class Grid {
   void get_dipole_grid_range(struct dipole_range &dr) const;
 
   // Helper function for set_interpolation_coefs
-  void set_interp_coef_sphere(const sphere_range &sr,
-                              const precision_t lon_in,
-                              const precision_t lat_in,
-                              const precision_t alt_in);
-  void set_interp_coef_cubesphere(const cubesphere_range &cr,
-                                  const precision_t lon_in,
-                                  const precision_t lat_in,
-                                  const precision_t alt_in);
+  struct interp_coef_t get_interp_coef_sphere(const sphere_range &sr,
+                                              const precision_t lon_in,
+                                              const precision_t lat_in,
+                                              const precision_t alt_in);
+  struct interp_coef_t get_interp_coef_cubesphere(const cubesphere_range &cr,
+                                                  const precision_t lon_in,
+                                                  const precision_t lat_in,
+                                                  const precision_t alt_in);
   // (note these are magnetic coordinates)
-  void set_interp_coef_dipole(const dipole_range &dr,
-                              const precision_t lon_in,
-                              const precision_t lat_in,
-                              const precision_t alt_in);
+  struct interp_coef_t get_interp_coef_dipole(const dipole_range &dr,
+					      const precision_t lon_in,
+					      const precision_t lat_in,
+					      const precision_t alt_in);
 
   // Processed interpolation coefficients
   std::vector<struct interp_coef_t> interp_coefs;

src/advance.cpp

@@ -46,6 +46,10 @@ bool advance(Planets &planet,
     didWork = neutralsMag.check_for_nonfinites("Top of Advance - ion grid");
   }
 
+  // here we are going to grab stuff from the neutral grid and put it on the
+  // ion grid
+  didWork = get_data_from_other_grid(gGrid, mGrid, neutrals.temperature_scgc, mGrid.test_scgc);
+
   json dummy = indices.get_all_indices(time.get_current());
 
   gGrid.calc_sza(planet, time);

src/grid.cpp

@@ -97,6 +97,9 @@ Grid::Grid(std::string gridtype) {
   geoAlt_scgc.set_size(nX, nY, nZ);
   geoLocalTime_scgc.set_size(nX, nY, nZ);
 
+  test_scgc.set_size(nX, nY, nZ);
+  test_scgc.zeros();
+
   refx_scgc.set_size(nX, nY, nZ);
   refy_scgc.set_size(nX, nY, nZ);
   refx_angle.set_size(nX, nY, nZ);