Merge branch 'develop' into cet

Author: abukowski21

doc/internals/coordinates.md

@@ -1,3 +1,35 @@
+
+# Coordinate Systems
+
+## GEO - Geographic coordinate system (or geodetic)
+- Longitude (radians) - radians east of prime meridian (0 - 2pi)
+- Latitude (radians) - angle between the equatorial plane and point. For oblate spheriod, this angle is not the same as the angle orthogonal to the surface of the planet and the equatorial plane.
+- Radius (meters) - the distance to the center of the planet (Altitude is often used instead).
+
+## PCPF - Planet-centered Planet-fixed or Geocentric coordinate system
+- Cartesian coordinates of Geographic coordinate system in meters (X, Y, Z).
+- X (meters) - aligned with the equator and prime meridian
+- Z (meters) - aligned with rotation axis of the planet
+- Y (meters) - completes the right-hand coordinate system
+
+## PSE - Planetary Solar Ecliptic coordinates
+- Cartesian coordinates tying together planet and Sun
+- X (meters) - points from the center of the planet to the Sun
+- Y (meters) - points from the center of the planet towards dusk and opposite planet's motion around the sun
+- Z (meters) - orthogonal to the ecliptic plane
+
+## Dipole Coordinates
+- Longitude (radians) - radians east of the meridian that contains the north magnetic pole and north rotation axis
+- P (?) - Identifies the field line, related to L-shell
+- Q (meters?) - The distance along the field line from some reference point, related to magnetic latitude.
+
+## More Dipole Coordinates
+- L-shell (Planetary Radii) - The distance from the planet's center at which the magnetic field encounters the dipole's equatorial plane
+- Magnetic Latitude (radians) - angle between the dipole's equatorial plane and the point.
+- Invariant Latitude (degrees) - angle between the dipole's equatorial plane and the point at which the field-line passes through a reference radius of the planet.  This is constant along the field-line and is related to the L-Shell.
+- Magnetic Local Time (hours) - Angle between the sun, the north magnetic pole, and the point. Explicitly, this is done in PSE XY coordinates, ignoring the Z coorinate.
+
+
 # Coordinates in Aether
 
 There are a variety of coordinates in Aether. This document describes some of them.
@@ -39,3 +71,4 @@ All locations should be described in the following coordinates:
 - lon, lat, radius (+alt)
 - magnetic lon, invariant lat?
 
+

include/calc_grid_derived.h

@@ -16,4 +16,12 @@ arma_vec calc_bin_edges(arma_vec centers);
 std::vector<precision_t> calc_bin_widths(std::vector<precision_t> centers);
 arma_vec calc_bin_widths(arma_vec centers);
 
+// ----------------------------------------------------------------------------
+// A helper function for mapping grids
+// ----------------------------------------------------------------------------
+bool grid_match(Grid gGrid, 
+                Grid mGrid, 
+                Quadtree gQuadtree,
+                Quadtree mQuadtree);
+
 #endif  // INCLUDE_CALC_GRID_DERIVED_H_

include/dipole.h

@@ -59,23 +59,23 @@ namespace Dipole2 {
  *************************************************/
 namespace Dipole4 {
 
-  /// The normalized origins of each face of the cube (i.e. corner)
+  /// The normalized origins of each node (i.e. corner)
   static const arma_mat ORIGINS = {
 				   { 0.0, -0.5, 0.0},
 				   { 1.0, -0.5, 0.0},
 				   { 1.0,  0.0, 0.0},
            { 0.0,  0.0, 0.0}
   };
 
-  /// Normalized right steps in cube
+  /// Normalized right steps in node
   static const arma_mat RIGHTS = {
 				  {1.0, 0.0, 0.0},
 				  {1.0, 0.0, 0.0},
 				  {1.0, 0.0, 0.0},
 				  {1.0, 0.0, 0.0}
   };
 
-  /// Normalized up steps in cube
+  /// Normalized up steps in node
   static const arma_mat UPS = {
 			       {0.0, 0.5, 0.0},
 			       {0.0, 0.5, 0.0},

include/grid.h

@@ -69,8 +69,15 @@ class Grid
   // These define the magnetic grid:
   // Armidillo Cube Versions:
   arma_cube magLon_scgc, magX_scgc;
+  // The magnetic latitude and altitude need to be defined better. This should be the angle between
+  // magnetic equator and the point, but sometimes it is invariant latitude.
   arma_cube magLat_scgc, magY_scgc;
+  // This is often just the altitude....
   arma_cube magAlt_scgc, magZ_scgc;
+  // Invariant latitude is the magnetic latitude that the field line hits at the lowest altitude.
+  // This is basically the L-shell, but models want it expressed as latitude and not L-shell.
+  arma_cube magInvLat_scgc;
+  // This is the angle from the sun, to the magnetic pole to the point.
   arma_cube magLocalTime_scgc;
 
   // Dipole coordinates:
@@ -152,7 +159,55 @@ class Grid
 
   // dr is the radial change along the third dimension, which is
   // primarily needed for building a hydrostatic solution
-  arma_cube dr_lower_scgc;
+  arma_cube dr_edge;
+
+  // i, j, k are the three directions, so these are the grid spacing
+  // between the cell centers in each direction, aligned with the grid
+  arma_cube i_center_scgc;
+  arma_cube j_center_scgc;
+  arma_cube k_center_scgc;
+
+  // edges are defined in the direction of the coordinate, shifted -half
+  // a cell in that direction.
+  arma_cube i_edge_scgc;
+  arma_cube j_edge_scgc;
+  arma_cube k_edge_scgc;
+
+  // corners are defined as shifted by -half a cell in each direction
+  arma_cube i_corner_scgc;
+  arma_cube j_corner_scgc;
+  arma_cube k_corner_scgc;
+
+  // native distances in native units:
+  arma_cube di_center_scgc;
+  arma_cube dj_center_scgc;
+  arma_cube dk_center_scgc;
+
+  // native distance in meters
+  arma_cube di_center_m_scgc;
+  arma_cube dj_center_m_scgc;
+  arma_cube dk_center_m_scgc;
+
+  // Gradients on the edges really only have to be between cells, so they
+  // can be defined at the interfaces (n-1 of them)
+  arma_cube di_edge;
+  arma_cube dj_edge;
+  arma_cube dk_edge;
+  // in meters:
+  arma_cube di_edge_m;
+  arma_cube dj_edge_m;
+  arma_cube dk_edge_m;
+
+  // These are for stretched grids:
+  arma_cube di_ratio;
+  arma_cube di_ratio_sq;
+  arma_cube di_one_minus_r2;
+  arma_cube dj_ratio;
+  arma_cube dj_ratio_sq;
+  arma_cube dj_one_minus_r2;
+  arma_cube dk_ratio;
+  arma_cube dk_ratio_sq;
+  arma_cube dk_one_minus_r2;
 
   arma_cube MeshCoefm2;
   arma_cube MeshCoefm1;
@@ -201,6 +256,7 @@ class Grid
   void set_variable_sizes();
 
   bool get_IsGeoGrid();
+  std::string get_gridtype();
   bool get_HasBField();
   void set_IsGeoGrid(bool value);
   void set_IsExperimental(bool value);
@@ -236,16 +292,21 @@ class Grid
   bool get_Is1Dy();
   bool get_Is1Dz();
 
-  void fill_grid(Planets planet);
+  //void fill_grid(Planets planet);
   void correct_xy_grid(Planets planet);
   void calc_sza(Planets planet, Times time);
   void calc_gse(Planets planet, Times time);
   void calc_mlt();
+  void calc_xyz(Planets planet);
 
   void calc_grid_spacing(Planets planet);
   void calc_alt_grid_spacing();
   void calc_lat_grid_spacing();
   void calc_long_grid_spacing();
+  void calc_maglong_grid_spacing();
+  void calc_i_grid_spacing();
+  void calc_j_grid_spacing();
+  void calc_k_grid_spacing();
 
   void fill_grid_radius(Planets planet);
   void calc_rad_unit(Planets planet);
@@ -389,6 +450,7 @@ class Grid
   bool IsExperimental;
   bool IsMagGrid;
   bool IsDipole = false;
+  std::string gridType;
 
   int64_t nX, nLons;
   int64_t nY, nLats;

include/ions.h

@@ -435,7 +435,7 @@ class Ions {
      \param dir the directory to read or write from/to
      \param DoRead whether to read (true) or write (false)
    **/
-  bool restart_file(std::string dir, bool DoRead);
+  bool restart_file(std::string dir, std::string cGridtype, bool DoRead);
 
   /**********************************************************************
      \brief Exchange messages between processors

include/neutrals.h

@@ -470,7 +470,7 @@ class Neutrals {
      \param dir directory to write restart files
      \param DoRead read the restart files if true, write if false
    **/
-  bool restart_file(std::string dir, bool DoRead);  
+  bool restart_file(std::string dir, std::string cGridtype, bool DoRead);  
 
   /**********************************************************************
      \brief Exchange messages between processors

include/solvers.h

@@ -91,6 +91,16 @@ std::vector<arma_cube> calc_gradient_cubesphere(arma_cube value, Grid grid);
 arma_cube calc_gradient_alt_4th(arma_cube value, Grid grid);
 arma_mat project_onesided_alt_3rd(arma_cube value, Grid grid, int64_t iAlt);
 
+// Calculate 4th-order gradients in the native coordinate system:
+arma_cube calc_gradient4o_i(arma_cube value, Grid grid);
+arma_cube calc_gradient4o_j(arma_cube value, Grid grid);
+arma_cube calc_gradient4o_k(arma_cube value, Grid grid);
+
+// Calculate 2nd-order gradients in the native coordinate system:
+arma_cube calc_gradient2o_i(arma_cube value, Grid grid);
+arma_cube calc_gradient2o_j(arma_cube value, Grid grid);
+arma_cube calc_gradient2o_k(arma_cube value, Grid grid);
+
 // interpolation in 1D
 precision_t linear_interpolation(const precision_t y0,
                                  const precision_t y1,

include/tools.h

@@ -65,6 +65,13 @@ std::string add_cmember(std::string inString);
 
 void display_vector(arma_vec vec);
 
+// ----------------------------------------------------------------------
+// Display an armadillo vector with a strong name in front
+// ----------------------------------------------------------------------
+
+void display_vector(std::string, arma_vec vec);
+void display_vector(std::string, std::vector<precision_t> vec);
+
 // ----------------------------------------------------------------------
 // synchronize a (boolean) variable across all processors
 // ----------------------------------------------------------------------

share/run/UA/inputs/defaults.json

@@ -100,8 +100,8 @@
 
     "Sources" : {
 	"Grid" : {
-	    "Centripetal" : true,
-		"Coriolis" : true },
+		"Coriolis" : true,
+		"Cent_acc": true },
 	"Neutrals" : {
 	    "NO_cool" : false,
 	    "O_cool": false },

src/advance.cpp

@@ -213,8 +213,12 @@ bool advance(Planets &planet,
 
     if (time.check_time_gate(input.get_dt_write_restarts())) {
       report.print(3, "Writing restart files");
-      neutrals.restart_file(input.get_restartout_dir(), DoWrite);
-      ions.restart_file(input.get_restartout_dir(), DoWrite);
+      neutrals.restart_file(input.get_restartout_dir(), gGrid.get_gridtype(),
+                            DoWrite);
+      neutralsMag.restart_file(input.get_restartout_dir(), mGrid.get_gridtype(),
+                               DoWrite);
+      ions.restart_file(input.get_restartout_dir(), gGrid.get_gridtype(), DoWrite);
+      ionsMag.restart_file(input.get_restartout_dir(), mGrid.get_gridtype(), DoWrite);
       time.restart_file(input.get_restartout_dir(), DoWrite);
     }
   }