Merge pull request #165 from AetherModel/ion_vertical_advection

Looks good! I added one more change - removing `centripetal` from the `.json` files since we are using `cent_acc` instead.

Author: abukowski21

doc/internals/coordinates.md

@@ -0,0 +1,74 @@
+
+# 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.
+
+## Spherical Coordinates
+
+The easiest coordinate system to understand within Aether is the spherical system, which is a longitude, latitude, radial (LLR) coordinate system. When the planet is a pure sphere, the LLR system is orthogonal - meaning that the grid lines up perfectly with the lines of constant longitude, latitude, and radius.
+
+In Aether, longitude and latitude are expressed in radians and are positive towards the east and towards the north. Radius is expressed in meters and is positive away from the planet (upwards).  In Aether, often Altitude is used instead of radius. When the planet is a perfect sphere, these are offset by a constant value.
+
+If the planet is an oblate spheriod, then the equator is larger than the pole, so that the planetary radius is dependent on latitude. Aether is currently set up so that Altitude is not dependent on latitude or longitude, so that if an oblate spheriod is used, then a constant altitude would have a radius that is dependent on latitude.  This means that the coordinate system is not purely orthogonal. At this time, this is not dealt with properly.
+
+Because Aether considers gravity to be a function of radius and explicitly includes the centrifugal acceleration, the pertubation away from a perfect sphere should mostly cancel.
+
+## i, j, k Coordinates
+
+As described in the grid.md file, Aether uses a logical '(i, j, k)' 3D grid structure.  Therefore, we refer to the 'primary' coordinates as the ijk coordinate system.  What this means is that the i-coordinate is in the i-direction, the j-coordinate is in the j-direction, and the k-coordinate is in the k-direction.  
+
+For the (perfectly) spherical grid, the i-coordinate is longitude, the j-coordinate is latitude, and the k-coordinate is radius.
+
+For the Cubedsphere grid, the i-coordinate is RIGHT, the j-coodinate is UP, and the k-coordinate is radius.  Each face of the cubedsphere has the same coordinate system, but only with reference to that face.  This means that if each face is looked at independently, the lower left corner is at (about) i = -45, j = -45 deg, while the upper right corner is at i = +45, j = +45 deg. Radius is treated the same as in a spherical grid.
+
+Should the official coordinates be in the native coordinates (which could be different for each system), or should the coordinates be in meters, such that when gradients are taken, they are in /m?
+
+Maybe we could have:
+
+i_scgc, j_scgc, k_scgc - coordinates in the native coordinates (radians, meters, etc.)
+im_scgc, jm_scgc, km_scgc - coordinate in meters
+
+The question is what variables do we need?
+
+Locations:
+- Cell Centers (these are the center of each volume)
+- Cell Edges in the i, j, k directions (these are the center of each area)
+- Cell Corners
+
+All locations should be described in the following coordinates:
+- i, j, k
+- lon, lat, radius (+alt)
+- magnetic lon, invariant lat?
+
+

ext/IE/EIE_IoLibrary.f90

@@ -281,6 +281,7 @@ subroutine IO_GetPotential(PotentialOut, iError)
   real, dimension(IOi_NeednMLTs,IOi_NeednLats)               :: ValueOut
   real, dimension(IOi_NeednMLTs,IOi_NeednLats), intent(out)  :: PotentialOut
   real :: Filler = 0.0
+  integer :: iLat, iMlt
 
   iError = 0
 
@@ -302,6 +303,27 @@ subroutine IO_GetPotential(PotentialOut, iError)
 
      call IO_GetNonGridBasedPotential(ValueOut, iError)
 
+    ! In Weimer, there are sometimes 0 potentials.  I don't know why.
+     if (maxval(abs(IOr2_NeedLats)) > 80.0) then
+
+      do iMlt = 2, IOi_NeednMLTs - 1
+        do iLat = 2, IOi_NeednLats - 1
+          if (ValueOut(iMlt, iLat) == 0.0) then
+            if (ValueOut(iMlt, iLat - 1)*ValueOut(iMlt, iLat + 1) > 0) then
+              ! this is a "hole" in the potential, since the potential is
+              ! the same sign both above and below the current point.
+              ValueOut(iMlt, iLat) = ( &
+                                     ValueOut(iMlt - 1, iLat) + &
+                                     ValueOut(iMlt, iLat - 1) + &
+                                     ValueOut(iMlt + 1, iLat) + &
+                                     ValueOut(iMlt, iLat + 1))/4.0
+            endif
+          endif
+        enddo
+      enddo
+
+     endif
+
      if (iError == 0) then
         PotentialOut = ValueOut
      else

include/advance.h

@@ -44,6 +44,7 @@ bool advance(Planets &planet,
 	     Chemistry &chemistry,
 	     Chemistry &chemistryMag,
 	     Electrodynamics &electrodynamics,
+	     Electrodynamics &electrodynamicsMag,
 	     Indices &indices,
 	     Logfile &logfile,
 	     Logfile &logfileMag);

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:
@@ -143,11 +150,65 @@ class Grid
   arma_cube sza_scgc;
   arma_cube cos_sza_scgc;
 
+  // dalt should be the altitudinal change along the third dimension,
+  // but is really the distance between grid points.
   arma_cube dalt_center_scgc;
   arma_cube dalt_lower_scgc;
   arma_cube dalt_ratio_scgc;
   arma_cube dalt_ratio_sq_scgc;
 
+  // dr is the radial change along the third dimension, which is
+  // primarily needed for building a hydrostatic solution
+  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;
   arma_cube MeshCoefp0;
@@ -195,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);
@@ -230,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);
@@ -383,6 +450,7 @@ class Grid
   bool IsExperimental;
   bool IsMagGrid;
   bool IsDipole = false;
+  std::string gridType;
 
   int64_t nX, nLons;
   int64_t nY, nLats;

include/inputs.h

@@ -436,13 +436,7 @@ class Inputs {
   bool get_O_cooling();
 
   /**********************************************************************
-     \brief returns settings["
-     \param 
-   **/
-  bool get_use_centripetal();
-  
-  /**********************************************************************
-     \brief returns settings["
+      \brief returns settings["
      \param 
    **/
   bool get_use_coriolis();
@@ -487,7 +481,8 @@ class Inputs {
   bool get_advection_neutrals_bulkwinds();
   bool get_advection_neutrals_implicitfriction();
 
-  
+  std::string get_advection_ions_along();
+
   /**********************************************************************
      \brief returns settings["
      \param